collegeSummary
Yakovlev D.M.
2024-12-20
options(warn = -1)
options(repr.plot.width=15, repr.plot.height=15)1 Предварительный анализ данных
library(readxl)
col_I_sn<-read_excel("stat2021_sm\\us_col.std\\I_shortname.xls")## New names:
## • `` -> `...1`col_I<-read_excel("stat2021_sm\\us_col.std\\I.xls")## New names:
## • `` -> `...1`head(col_I_sn)|>rmarkdown::paged_table()summary(col_I_sn)## ...1 PPIND FICE STATE
## Length:176 Min. :1.000 Min. : 1009 Length:176
## Class :character 1st Qu.:1.000 1st Qu.: 1738 Class :character
## Mode :character Median :1.000 Median : 2550 Mode :character
## Mean :1.301 Mean : 2901
## 3rd Qu.:2.000 3rd Qu.: 3417
## Max. :2.000 Max. :10366
##
## TYPE AVRMATH AVRVERB AVRCOMB
## Length:176 Min. :390.0 Min. :391.0 Min. : 810
## Class :character 1st Qu.:520.0 1st Qu.:454.5 1st Qu.: 980
## Mode :character Median :544.0 Median :479.5 Median :1017
## Mean :563.5 Mean :494.7 Mean :1058
## 3rd Qu.:603.0 3rd Qu.:518.5 3rd Qu.:1126
## Max. :750.0 Max. :665.0 Max. :1410
## NA's :68 NA's :68 NA's :67
## AVR_ACT MATH_1 MATH_3 VERB_1
## Min. :19.00 Min. :350.0 Min. :480.0 Min. :320.0
## 1st Qu.:22.00 1st Qu.:460.0 1st Qu.:590.0 1st Qu.:400.0
## Median :23.00 Median :500.0 Median :627.5 Median :430.0
## Mean :23.69 Mean :515.4 Mean :635.0 Mean :448.3
## 3rd Qu.:25.00 3rd Qu.:570.0 3rd Qu.:687.5 3rd Qu.:480.0
## Max. :31.00 Max. :740.0 Max. :780.0 Max. :630.0
## NA's :83 NA's :34 NA's :34 NA's :34
## VERB_3 ACT_1 ACT_3 APP_REC
## Min. :440.0 Min. :16.00 Min. :21.00 Min. : 787
## 1st Qu.:520.0 1st Qu.:19.00 1st Qu.:25.00 1st Qu.: 4323
## Median :550.0 Median :21.00 Median :26.00 Median : 7654
## Mean :561.3 Mean :21.28 Mean :26.55 Mean : 8544
## 3rd Qu.:600.0 3rd Qu.:23.00 3rd Qu.:28.00 3rd Qu.:11776
## Max. :720.0 Max. :29.00 Max. :33.00 Max. :48094
## NA's :34 NA's :67 NA's :67 NA's :1
## APP_ACC NEW_STUD NEW10 NEW25 FULLTIME
## Min. : 507 Min. : 210 Min. : 8.00 Min. :24.00 Min. : 912
## 1st Qu.: 3033 1st Qu.:1264 1st Qu.:24.00 1st Qu.:52.00 1st Qu.: 5846
## Median : 4761 Median :1949 Median :32.00 Median :63.00 Median :10215
## Mean : 5546 Mean :2252 Mean :41.48 Mean :65.55 Mean :11296
## 3rd Qu.: 7232 3rd Qu.:3035 3rd Qu.:57.00 3rd Qu.:82.00 3rd Qu.:15033
## Max. :26330 Max. :7425 Max. :98.00 Max. :99.00 Max. :31643
## NA's :1 NA's :16 NA's :26
## PARTTIME IN_STATE OUT_STAT R_B_COST
## Min. : 16.0 Min. : 647 Min. : 2279 Min. :2082
## 1st Qu.: 804.8 1st Qu.: 2100 1st Qu.: 6712 1st Qu.:3588
## Median : 1694.5 Median : 3030 Median : 8400 Median :4213
## Mean : 2519.9 Mean : 6641 Mean :10108 Mean :4538
## 3rd Qu.: 3240.8 3rd Qu.:12348 3rd Qu.:12668 3rd Qu.:5564
## Max. :21836.0 Max. :20100 Max. :20100 Max. :7425
## NA's :6 NA's :7 NA's :1
## ROOM BOARD ADD_FEE BOOK PERSONAL
## Min. :1033 Min. :1000 Min. : 20.0 Min. : 300.0 Min. : 300
## 1st Qu.:1810 1st Qu.:1762 1st Qu.: 210.0 1st Qu.: 500.0 1st Qu.:1164
## Median :2644 Median :2125 Median : 425.5 Median : 600.0 Median :1600
## Mean :2708 Mean :2241 Mean : 648.1 Mean : 603.1 Mean :1763
## 3rd Qu.:3490 3rd Qu.:2586 3rd Qu.: 694.0 3rd Qu.: 673.8 3rd Qu.:2200
## Max. :6965 Max. :4760 Max. :4374.0 Max. :1230.0 Max. :6800
## NA's :42 NA's :64 NA's :40 NA's :2 NA's :11
## PH_D TERM_D SF_RATIO DONATE
## Min. :63.00 Min. :67.00 Min. : 2.90 Min. : 4.00
## 1st Qu.:80.50 1st Qu.:87.00 1st Qu.:10.88 1st Qu.:10.75
## Median :87.00 Median :92.00 Median :14.50 Median :17.00
## Mean :85.72 Mean :90.52 Mean :14.23 Mean :19.01
## 3rd Qu.:92.00 3rd Qu.:96.00 3rd Qu.:18.02 3rd Qu.:24.00
## Max. :99.00 Max. :99.00 Max. :24.70 Max. :54.00
## NA's :9 NA's :16 NA's :12
## INSTRUCT GRADUAT SAL_FULL SAL_AC
## Min. : 3605 Min. :10.00 Min. : 446.0 Min. :364.0
## 1st Qu.: 7604 1st Qu.:47.50 1st Qu.: 597.0 1st Qu.:445.0
## Median : 9840 Median :62.00 Median : 661.0 Median :479.5
## Mean :12832 Mean :62.02 Mean : 669.2 Mean :487.1
## 3rd Qu.:14340 3rd Qu.:74.50 3rd Qu.: 732.2 3rd Qu.:521.2
## Max. :62469 Max. :99.00 Max. :1009.0 Max. :733.0
## NA's :5
## SAL_AS SAL_ALL COMP_FUL COMP_AC
## Min. :323.0 Min. :362.0 Min. : 537.0 Min. :438.0
## 1st Qu.:381.8 1st Qu.:472.2 1st Qu.: 729.0 1st Qu.:556.5
## Median :407.0 Median :522.5 Median : 815.0 Median :606.0
## Mean :412.8 Mean :534.0 Mean : 827.2 Mean :612.1
## 3rd Qu.:434.2 3rd Qu.:578.2 3rd Qu.: 910.0 3rd Qu.:662.2
## Max. :576.0 Max. :866.0 Max. :1236.0 Max. :909.0
##
## COMP_AS COMP_ALL NUM_FULL NUM_AC
## Min. :395.0 Min. : 436.0 Min. : 39.0 Min. : 32.0
## 1st Qu.:481.0 1st Qu.: 587.0 1st Qu.:184.5 1st Qu.:138.5
## Median :508.5 Median : 652.0 Median :278.0 Median :208.0
## Mean :519.0 Mean : 665.8 Mean :336.8 Mean :230.3
## 3rd Qu.:552.2 3rd Qu.: 730.2 3rd Qu.:457.8 3rd Qu.:299.0
## Max. :717.0 Max. :1075.0 Max. :997.0 Max. :721.0
## NA's :1
## NUM_AS NUM_INS NUM_ALL
## Min. : 29.0 Min. : 0.00 Min. : 108.0
## 1st Qu.:128.5 1st Qu.: 5.00 1st Qu.: 505.5
## Median :175.0 Median : 16.00 Median : 721.0
## Mean :190.7 Mean : 27.09 Mean : 812.8
## 3rd Qu.:238.2 3rd Qu.: 40.00 3rd Qu.:1035.0
## Max. :510.0 Max. :178.00 Max. :2261.0
## NA's :1col_I_sn$PPIND=ifelse(col_I_sn$PPIND==1, "public", "private")NA review
col_I_sn_noNA<-col_I_sn|>na.omit()
paste("col_I_sn dim: ", paste(dim(col_I_sn),collapse=", "))## [1] "col_I_sn dim: 176, 49"paste("col_I_sn_noNA dim: ", paste(dim(col_I_sn_noNA), collapse=", "))## [1] "col_I_sn_noNA dim: 22, 49"naCol<-colSums(is.na(col_I_sn))|>sort(decreasing=TRUE)
naCol## AVR_ACT AVRMATH AVRVERB AVRCOMB ACT_1 ACT_3 BOARD ROOM
## 83 68 68 67 67 67 64 42
## ADD_FEE MATH_1 MATH_3 VERB_1 VERB_3 NEW25 NEW10 TERM_D
## 40 34 34 34 34 26 16 16
## DONATE PERSONAL PH_D IN_STATE PARTTIME GRADUAT BOOK APP_REC
## 12 11 9 7 6 5 2 1
## APP_ACC OUT_STAT NUM_AC NUM_INS ...1 PPIND FICE STATE
## 1 1 1 1 0 0 0 0
## TYPE NEW_STUD FULLTIME R_B_COST SF_RATIO INSTRUCT SAL_FULL SAL_AC
## 0 0 0 0 0 0 0 0
## SAL_AS SAL_ALL COMP_FUL COMP_AC COMP_AS COMP_ALL NUM_FULL NUM_AS
## 0 0 0 0 0 0 0 0
## NUM_ALL
## 0table(naCol)## naCol
## 0 1 2 5 6 7 9 11 12 16 26 34 40 42 64 67 68 83
## 21 5 1 1 1 1 1 1 1 2 1 4 1 1 1 3 2 1log dataset
q_cols<-sapply(col_I_sn, is.numeric)
q_cols[c("PPIND", "FICE")]<-FALSE
q_cols<-names(q_cols)[q_cols==TRUE]
col_I_sn_log<-col_I_sn|>
mutate(across(q_cols, ~log(.x)))|>
rename_with(~paste0("log-", .x), q_cols)
head(col_I_sn_log)|>rmarkdown::paged_table()1.1 Разобраться в том, что означают признаки.
Выделим основные признаки, по которым будет проходить дальнейший анализ данных.
1.1.1 Список признаков.
1.1.1.1 Количественные признаки:
- AVRMATH Average Math SAT score
- AVRVERB Average Verbal SAT score
- AVRCOMB Average Combined SAT score
- AVR_ACT Average ACT score
- MATH_1 First quartile - Math SAT
- MATH_3 Third quartile - Math SAT
- VERB_1 First quartile - Verbal SAT
- VERB_3 Third quartile - Verbal SAT
- ACT_1 First quartile - ACT
- ACT_3 Third quartile - ACT
- APP_REC Number of applications received
- APP_ACC Number of applicants accepted
- NEW_STUD Number of new students enrolled
- FULLTIME Number of fulltime undergraduates
- PARTTIME Number of parttime undergraduates
- IN_STATE In-state tuition
- OUT_STAT Out-of-state tuition
- R_B_COST Room and board costs
- ROOM Room costs
- BOARD Board costs
- ADD_FEE Additional fees
- BOOK Estimated book costs
- PERSONAL Estimated personal spending
- PH_D Pct. of faculty with Ph.D.’s
- TERM_D Pct. of faculty with terminal degree
- SAL_FULL Average salary - full professor
- SAL_AC Average salary - associate professor
- SAL_AS Average salary - assistant professor
- SAL_ALL Average salary - all ranks
- COMP_FUL Average compensation - full professor
- COMP_AC Average compensation - associate professor
- COMP_AS Average compensation - assistant professor
- COMP_ALL Average compensation - all ranks
- NUM_FULL Number of full professor
- NUM_AC Number of associate professor
- NUM_AS Number of assistant professor
- NUM_INS Number of instructors
- NUM_ALL Number of faculty - all ranks
- INSTRUCT Instructional expenditure per student
- GRADUAT Graduation rate
- SF_RATIO Student/faculty ratio
- DONATE Pct.alumni who donate
- NEW10 Pct. new students from top 10% of H.S. class - % студентов из топ 10% своей старшей школы
- NEW25 Pct. new students from top 25% of H.S. class - % студентов из топ 25% своей старшей школы
1.1.1.2 Качественные признаки:
- FICE - Federal ID Number
- …1 - Название университета
- PPIND Public/private indicator (public=1, private=2)
- STATE State (postal code)
- TYPE - I (можно удалить)
Всего 49 признаков.
1.2 Если признаков очень много, то отобрать признаки (примерно 7-10)…
Отберём из условий задачи и по смыслу
- PPIND (для группировки)
Из условий задачи: 2. ADD_FEE - дополнительные выплаты 3. BOOK - плата за бронирование 4. NEW10 (зависимая переменная, также NEW25)
NEW10 зависит от AVRMATH, MATH_1, MATH_3, AVRVERB, VERB_1, VERB_3, AVR_ACT, ACT_1, ACT_3 и AVRCOMB. Признаков много, так что выберем самый обобщающий - AVRCOMB 5. AVRCOMB
SF_RATIO
PH_D
GRADUAT 6-8 из предположения, что они взаимосвязаны
R_B_COST
IN_STATE
OUT_STAT 9-11 - зависимость стоимости проживания (как окажется дальше из matrix plot, IN_STATE - фиктивный признак, то есть не информативный)
1.3 Определить вид признаков (колич., порядковые, качеств.)…
Проделаем это для рассматриваемых 11 признаков.
Setup
colNames<-c("ADD_FEE", "BOOK", "R_B_COST","IN_STATE","OUT_STAT","NEW10","AVRCOMB","SF_RATIO","PH_D","GRADUAT","PPIND")
colFin<-c("ADD_FEE", "BOOK", "PPIND")
colCost<-c("R_B_COST", "IN_STATE", "OUT_STAT", "PPIND")
colNew<-c("NEW10", "AVRCOMB", "PPIND")
colStud<-c("SF_RATIO", "PH_D", "GRADUAT", "PPIND")
col_I_comp<-col_I_sn[colNames]
col_I_comp_log<-col_I_sn_log[c(paste0("log-", colNames)[-length(colNames)], "PPIND")]
head(col_I_comp_log)|>rmarkdown::paged_table()unique_info <- function(column) {
column<-na.omit(column)
list(ratio=length(unique(column)) / length(column), unique=length(unique(column)), total=length(column), moda = max(table(column)))
}
ratios <- sapply(col_I_comp, unique_info)
ratios## ADD_FEE BOOK R_B_COST IN_STATE OUT_STAT NEW10 AVRCOMB
## ratio 0.875 0.3793103 0.9431818 0.9053254 0.9314286 0.425 0.7889908
## unique 119 66 166 153 163 68 86
## total 136 174 176 169 175 160 109
## moda 4 29 3 4 4 10 3
## SF_RATIO PH_D GRADUAT PPIND
## ratio 0.6363636 0.1976048 0.3918129 0.01136364
## unique 112 33 67 2
## total 176 167 171 176
## moda 5 10 7 1231.3.1 Количественные признаки
1.3.1.1 Непрерывные признаки
- ADD_FEE
- R_B_COST
- IN_STATE
- OUT_STATE
- AVRCOMB
- SF_RATIO
- GRADUAT
1.3.1.2 Дискретные признаки
- BOOK (непрерывный, но дискретный из-за моды. Может быть округление)
- PH_D (считаем из-за плохой точности и округления данных)
- NEW10 (по ratio и moda, низкая точность)
1.3.2 Качественные признаки
- PPIND - private/public.
1.4 Построить matrix plot (pairs plot), его долго разглядывать с точки зрения outliers, неоднородностей, вида распределений, вида зависимостей (линейные/нелинейные) и пр.
1.4.1 Первое построение pairs plot. Общие наблюдения.
ggpairs(col_I_comp, aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))
1.4.1.1 Наблюдения:
- IN_STATE и OUT_STAT для private выстраивает линейную зависимость, для public - близкая к линейной зависимости. Отличительная особенность признака IN_STATE от OUT_STAT - наличие других факторов на образование цены обучения. OUT_STATE более приложим для анализа данных, поэтому IN_STATE можно считать фиктивным признаком.
- ADD_FEE похоже на log-normal. Заметны Outliers.
1.4.1.1.1 Outliers, ADD_FEE:
col_I_sn|>filter(ADD_FEE>3000)University of California, Massachusetts - особые индивиды, в связи с высоким ВВП штатов.
- По histogram plot неоднородности (по PPIND) в R_B_COST, IN_STATE, OUT_STAT, SF_RATIO, GRADUAT.
Видим, что NEW10 зависит от PH_D, GRADUAT, AVRCOMB, отдельно для private зависит от OUT_STAT, IN_STATE, для public - от ADD_FEE, BOOK, R_B_COST
1.5 Если есть сильно несимметричные (с хвостом вправо) распределения на положительной полуоси, то…
1.5.1 Второе построение pairs plot.
Прологарифмируем ADD_FEE
col_I_mixed<-col_I_comp|>mutate(ADD_FEE=col_I_comp_log$`log-ADD_FEE`)|>rename(`log-ADD_FEE`=ADD_FEE)|>dplyr::select(-IN_STATE)
head(col_I_mixed)ggpairs(col_I_mixed, aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
1.5.1.1 Наблюдения:
- Неоднородности (по PPIND) R_B_COST, OUT_STAT, NEW10.
- Близкие к линейномй зависимости NEW10/ AVRCOMB.
- Outliers, BOOK
1.5.1.1.1 Outliers, BOOK
Особые индивиды
col_I_sn|>filter(BOOK>=900 | BOOK <= 350)Howard, Tulane, Hofstra University, Rensselaer Polytechn - особенности расположения, штата.
1.6 Если есть outliers, то попробовать объяснить причину (ошибка в данных, особые индивиды) и удалить их.
1.6.1 BOOK
ggpairs(col_I_mixed[,c("log-ADD_FEE", "BOOK", "PPIND")], aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
col_I_sn|>filter(BOOK>=900 | BOOK <= 350)Может быть как низкой, так и высокой в зависимости от штата, политики университета и скрытых факторов.
1.6.2 ADD_FEE
col_I_sn|>filter(ADD_FEE>3000)University of California, Massachusetts - особые индивиды (из-за ВВП).
1.6.3 AVRCOMB/NEW10
ggpairs(col_I_mixed[colNew], aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
col_I_sn|>filter(AVRCOMB<=1200 & NEW10>=75)University of California, North Carolina. Причиной может быть низкий порог для поступления.
1.6.4 PH_D
ggpairs(col_I_mixed[colStud], aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
col_I_sn|>filter(PH_D<=70)Незначительные выбросы - особые индивиды. Могут быть связаны с низким GRADUAT и высоким SF_RATIO.
1.6.5 OUT_STAT
ggpairs(col_I_mixed[,c("R_B_COST", "OUT_STAT", "PPIND")], aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
col_I_sn|>filter(OUT_STAT>=14000 & PPIND=="public")col_I_sn|>filter(OUT_STAT<=10500 & PPIND=="private")Цены на обучение варьируются по штатам и самим университетам.
1.7 Если есть неоднородности (например, видны два облака точек), то объяснить причину (найти категоризующую переменную, объясняющую эту неоднородность).
ggpairs(col_I_mixed, aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Выделим неоднородности R_B_COST, OUT_STAT. Явное разделение на два облака точек для этих признаков может быть связано с разницей финансирования. Для public - государство, штат и граждане, для private - только граждане.
1.8 В дальнейшем вид pairs plots, распределения признаков и корреляции анализировать отдельно для неоднородных групп.
1.9 Посмотрите также на descriptive statistics с точки зрения минимумов-максимумов, асимметрии, эксцесса и пр.
results <- describeBy(col_I_mixed, group=col_I_comp$PPIND)
df_results <- map_dfr(results, ~as.data.frame(.x), .id = "group")
kable(df_results, "html") %>%
kable_styling(full_width = FALSE) %>%
scroll_box(width = "100%", height = "400px")| group | vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| log-ADD_FEE…1 | private | 1 | 46 | 5.838785 | 0.7936613 | 6.003811 | 5.872662 | 0.7083206 | 3.688880 | 7.515345 | 3.826465 | -0.4493357 | 0.0057261 | 0.1170190 |
| BOOK…2 | private | 2 | 53 | 610.075472 | 158.1109456 | 600.000000 | 597.534884 | 148.2600000 | 300.000000 | 1230.000000 | 930.000000 | 1.2325851 | 3.0965764 | 21.7182087 |
| R_B_COST…3 | private | 3 | 53 | 5839.811321 | 997.0704793 | 5975.000000 | 5922.255814 | 816.9126000 | 3320.000000 | 7425.000000 | 4105.000000 | -0.7233444 | -0.2295138 | 136.9581633 |
| OUT_STAT…4 | private | 4 | 53 | 15747.358491 | 3924.9053349 | 17020.000000 | 16306.186047 | 2772.4620000 | 2340.000000 | 20100.000000 | 17760.000000 | -1.3082365 | 1.3192370 | 539.1272103 |
| NEW10…5 | private | 5 | 52 | 57.788461 | 24.3796255 | 58.000000 | 57.761905 | 33.3585000 | 16.000000 | 98.000000 | 82.000000 | 0.0381849 | -1.4295240 | 3.3808458 |
| AVRCOMB…6 | private | 6 | 29 | 1177.034483 | 144.2860459 | 1199.000000 | 1177.680000 | 161.6034000 | 883.000000 | 1410.000000 | 527.000000 | -0.0893581 | -1.2296542 | 26.7932461 |
| SF_RATIO…7 | private | 7 | 53 | 9.813208 | 4.3268854 | 9.200000 | 9.555814 | 4.8925800 | 2.900000 | 20.500000 | 17.600000 | 0.4758614 | -0.7282486 | 0.5943434 |
| PH_D…8 | private | 8 | 49 | 89.591837 | 7.8658502 | 91.000000 | 90.317073 | 7.4130000 | 71.000000 | 99.000000 | 28.000000 | -0.8078233 | -0.4898133 | 1.1236929 |
| GRADUAT…9 | private | 9 | 52 | 78.711539 | 15.7534142 | 78.500000 | 80.023809 | 17.7912000 | 33.000000 | 99.000000 | 66.000000 | -0.6127581 | -0.1790876 | 2.1846055 |
| PPIND…10 | private | 10 | 53 | 1.000000 | 0.0000000 | 1.000000 | 1.000000 | 0.0000000 | 1.000000 | 1.000000 | 0.000000 | NaN | NaN | 0.0000000 |
| log-ADD_FEE…11 | public | 1 | 90 | 6.014273 | 1.1296149 | 6.093567 | 6.024473 | 0.8723448 | 2.995732 | 8.383433 | 5.387701 | -0.1980167 | 0.2896088 | 0.1190719 |
| BOOK…12 | public | 2 | 121 | 600.008264 | 101.1621220 | 600.000000 | 591.103093 | 111.1950000 | 400.000000 | 858.000000 | 458.000000 | 0.6447669 | -0.2889320 | 9.1965565 |
| R_B_COST…13 | public | 3 | 123 | 3977.073171 | 842.4904744 | 3811.000000 | 3911.979798 | 722.0262000 | 2082.000000 | 6607.000000 | 4525.000000 | 0.7886999 | 0.7902446 | 75.9648078 |
| OUT_STAT…14 | public | 4 | 122 | 7657.549180 | 2252.3124723 | 7446.500000 | 7517.163265 | 1939.2408000 | 2279.000000 | 15732.000000 | 13453.000000 | 0.7771996 | 1.2996173 | 203.9147900 |
| NEW10…15 | public | 5 | 108 | 33.620370 | 21.8669373 | 26.000000 | 29.920455 | 11.8608000 | 8.000000 | 95.000000 | 87.000000 | 1.5420644 | 1.5799521 | 2.1041470 |
| AVRCOMB…16 | public | 6 | 80 | 1014.662500 | 84.6586510 | 997.000000 | 1010.687500 | 75.6126000 | 810.000000 | 1240.000000 | 430.000000 | 0.4180047 | 0.1997793 | 9.4651249 |
| SF_RATIO…17 | public | 7 | 123 | 16.129268 | 3.8302021 | 16.500000 | 16.178788 | 4.1512800 | 6.700000 | 24.700000 | 18.000000 | -0.0994206 | -0.5699091 | 0.3453577 |
| PH_D…18 | public | 8 | 118 | 84.110169 | 7.7985139 | 85.000000 | 84.406250 | 7.4130000 | 63.000000 | 99.000000 | 36.000000 | -0.3698537 | -0.3490528 | 0.7179114 |
| GRADUAT…19 | public | 9 | 119 | 54.731092 | 15.5730632 | 54.000000 | 54.412371 | 16.3086000 | 10.000000 | 95.000000 | 85.000000 | 0.1127780 | -0.1892703 | 1.4275804 |
| PPIND…20 | public | 10 | 123 | 2.000000 | 0.0000000 | 2.000000 | 2.000000 | 0.0000000 | 2.000000 | 2.000000 | 0.000000 | NaN | NaN | 0.0000000 |
1.9.1 Наблюдения:
1.9.1.1 Mean, median, skew:
NEW10 и AVRCOMB имеют skew, близкий к 0 (private), SF_RATIO, GRADUAT, log-ADD_FEE околонулевой skew (public).
1.9.1.2 Kurtosis:
Близкий к 0 (private): log-ADD_FEE (но skew ~ -0.44).
2 2. О виде распределений и о сравнении распределений
2.1 Первые два пункта индивидуального задания нужно делать не по указанному порядку, а как того требует логика статистики. Чтобы сравнивать выборки по t-критерию, нужно знать о том, близки ли распределения в сравниваемых группах к нормальным или хотя бы к унимодальным и симметричным. Чтобы проверять распределения признаков на нормальность, нужно знать, что рассматривается однородная выборка.
col_split_mixed<-split(col_I_mixed, col_I_mixed$PPIND)
names(col_split_mixed)<-c("private", "public")
colMixed<-names(col_I_mixed)
colMixed<-colMixed[-length(colMixed)]
colMixed## [1] "log-ADD_FEE" "BOOK" "R_B_COST" "OUT_STAT" "NEW10"
## [6] "AVRCOMB" "SF_RATIO" "PH_D" "GRADUAT"2.1.1 Проверка на нормальность по хи-квадрат критерию (Авторство: Яковлев Д.М.)
chi_squared_normality_test <- function(data) {
data <- na.omit(data)
n <- length(data)
if (n < 10) stop("Sample size too small for this test. Need at least 10 observations.")
mu <- mean(data)
sigma <- sd(data)
# Calculate quantiles for breaks
probs <- seq(5/n, 1 - 5/n, by = 5/n)
breaks <- qnorm(probs, mean = mu, sd = sigma)
breaks <- c(qnorm(0), breaks, qnorm(1))
#Observed Frequencies
observed <- hist(data, breaks = breaks, plot = FALSE)$counts
#Expected Frequencies (using pnorm)
expected <- diff(pnorm(breaks, mean = mu, sd = sigma)) * n
chi_squared <- sum((observed - expected)^2 / expected)
df <- length(expected) - 3 #Corrected degrees of freedom
p_value <- pchisq(chi_squared, df = df, lower.tail = FALSE)
return(list(p_value = p_value, statistic = chi_squared, df = df))
}2.1.1.1 Проверка на достоверность (p-value должно иметь равномерное распределение на (0, 1))
n = 500
test<-replicate(5000, chi_squared_normality_test(rnorm(n))$p_value)
quantile(test)## 0% 25% 50% 75% 100%
## 0.000116153 0.249925966 0.509619748 0.761571168 0.9996798142.1.2 Тесты на нормальность
perform_normality_tests <- function(df, cols) {
# Input validation (remains the same)
# Function to perform a single normality test (simplified)
run_test <- function(data, test_name) {
if (test_name == "Chi-squared") {
chisq_test<-chi_squared_normality_test(data)
return(chisq_test[c("p_value", "statistic", "df")])
}
test_result <- switch(test_name,
"Lilliefors" = lillie.test(data),
"Anderson-Darling" = ad.test(data),
"Shapiro-Wilk" = shapiro.test(data)
)
return(list(p_value = test_result$p.value, statistic = test_result$statistic, df=length(na.omit(data)))) # Removed test_name
}
results <- lapply(cols, function(col) {
data <- df[[col]]
n_na <- sum(is.na(data))
if (n_na > 0) {
warning(paste0("Removed ", n_na, " NA values from column '", col, "' before testing."))
data <- na.omit(data)
}
test_results <- lapply(c("Lilliefors", "Anderson-Darling", "Shapiro-Wilk", "Chi-squared"), function(test) run_test(data, test))
list(column = col, tests = test_results)
})
results_df <- do.call(rbind, lapply(results, function(x) {
test_names <- c("Lilliefors", "Anderson-Darling", "Shapiro-Wilk", "Chi-squared")
p_values <- sapply(x$tests, "[[", "p_value")
data.frame(
Column = x$column,
Test = test_names, #Directly use test_names vector
Statistic = sapply(x$tests, "[[", "statistic"),
P_value = p_values,
Size = sapply(x$tests, "[[", "df"),
Significance = sapply(p_values, function(p) ifelse(p < 0.05, "\u2714", "\u2716")),
stringsAsFactors = FALSE
)
}))
plots <- lapply(cols, function(col) {
data <- df[[col]]
data <- na.omit(data)
ggqqplot(data, title = paste("Normal Probability Plot for", col))
})
names(plots)<-cols
list(results = results_df, plots = plots)
}2.2 Так как визуально однородность при предварительном анализе была уже исследована, то можно начинать с анализа вида распределения признаков, возможно, по группам. Сюда входит: normal probability plot (что это такое?), проверка по критериям Лиллиефорса, AD, хи-квадрат, Шапиро-Уилка. По критерию хи-квадрат, а также визуально по PP-plot можно проверить и гипотезы о согласии с другими распределениями, например, логнормальным. Задаваемые вопросы: чем отличается критерий Лиллиефорса от критерия Колмогорова, в чем отличие AD, как выглядят статистики критериев, что такое PP-plot и normal probability plot, почему естественно при рисовании normal probability plot одновременно смотреть на результаты критерия Шапиро-Уилка.
2.2.1 Public: проверка на нормальность, тесты
# colSums(is.na(col_split_mixed$public))
col_norm_public<-col_split_mixed$public|>perform_normality_tests(colMixed)
col_norm_public$plots## $`log-ADD_FEE`
##
## $BOOK
##
## $R_B_COST
##
## $OUT_STAT
##
## $NEW10
##
## $AVRCOMB
##
## $SF_RATIO
##
## $PH_D
##
## $GRADUAT
col_norm_public$results|>dplyr::select(-Statistic)|>rmarkdown::paged_table()2.2.1.1 Наблюдения:
Хи-квадрат: не отвергаем AVRCOMB, SF_RATIO, PH_D, GRADUAT, R_B_COST, OUT_STAT,
Можно заметить “хвосты” на некоторых графиках.
2.2.2 Private: проверка на нормальность, тесты
col_norm_private<-col_split_mixed$private|>perform_normality_tests(colMixed)
col_norm_private$plots## $`log-ADD_FEE`
##
## $BOOK
##
## $R_B_COST
##
## $OUT_STAT
##
## $NEW10
##
## $AVRCOMB
##
## $SF_RATIO
##
## $PH_D
##
## $GRADUAT
col_norm_private$results|>dplyr::select(-Statistic)|>rmarkdown::paged_table()2.2.2.1 Наблюдения:
Хи-квадрат: не отвергаем log-ADD_FEE, R_B_COST, NEW10, AVRCOMB, SF_RATIO, GRADUAT
2.3 Сначала имеет смысл посмотреть на сравнение распределений в группах с помощью ящиков с усами. С помощью ящиков с усами там, где групп больше двух, можно выбрать две из них, которые интересно сравнить с помощью критериев.
Совместно со следующим заданием.
2.4 Если в задании есть сравнение независимых выборок (точнее, распределений независимых случайных величин), то начинать нужно с t-критерия, который мощный против альтернатив, заключающихся в наиболее легко интерпретируемом сдвиге (т.е. разнице средних). Нужно не забыть, что у критерия есть варианты для модели с одинаковыми дисперсиями (получается точное p-value, которое может быть неправильным, если на самом деле дисперсии одинаковые) и с произвольными дисперсиями. В результате, получатся два критерия для гипотезы о равенстве средних и два критерия о равенстве разбросов. Нужно уметь объяснять, что это за критерии и при каких условиях их можно применять. Не забудьте, что при использовании асимптотических критериев нужно обращать внимание на объемы выборок. Сделайте выводы о том, для каких признаков есть разница в сдвиге.
2.4.1 Тесты на равенство дисперсий/матожиданий, тест Wilcoxon.
perform_tests <- function(data,
quantitative_columns,
categorical_column,
alpha = 0.05) {
results <- list()
for (col in quantitative_columns) {
levene_test<-leveneTest(data[[col]] ~ data[[categorical_column]])
var_test <- var.test(data[[col]] ~ data[[categorical_column]])
wilcox_test <- wilcox.test(data[[col]]~data[[categorical_column]])
t_test_equal <- t.test(data[[col]] ~ data[[categorical_column]], var.equal = TRUE)
t_test_unequal <- t.test(data[[col]] ~ data[[categorical_column]], var.equal = FALSE)
ks_test<-ks.test(data[[col]]~data[[categorical_column]])
results[[col]] <- list(
t_test_unequal_p_value = t_test_unequal$p.value,
t_test_equal_p_value = t_test_equal$p.value,
wilcox_test_p_value = wilcox_test$p.value,
levene_test_p_value = levene_test$`Pr(>F)`[1],
var_test_p_value = var_test$p.value,
ks_test_p_value = ks_test$p.value
)
}
return(results |> enframe(name = "attribute", value = "p_values") |>
unnest_wider(p_values))
}df_tests = perform_tests(col_I_mixed, colMixed, "PPIND")df_long <- col_I_mixed|>
pivot_longer(cols = -PPIND, names_to = "attribute", values_to = "value")# Объединяем результаты теста с основными данными
df_with_tests <- df_long %>%
left_join(df_tests, by = "attribute")df_with_tests_cost <- df_with_tests[df_with_tests$attribute%in%colCost[-length(colCost)],]
df_with_tests_new <- df_with_tests[df_with_tests$attribute%in%colNew[-length(colNew)],]
df_with_tests_stud <- df_with_tests[df_with_tests$attribute%in%colStud[-length(colStud)],]
df_with_tests_fin <- df_with_tests[df_with_tests$attribute%in%c("log-ADD_FEE", "BOOK"),]
p<-ggplot(df_with_tests_cost,
aes(x = attribute,
y = value,
fill = PPIND)) +
geom_boxplot() +
labs(title = "Boxplot",
x = "Attribute",
y = "Value")+
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1), plot.caption = element_text(size = 8, hjust = 0)) +
scale_fill_manual(values = c("private" = "red", "public" = "lightblue"))2.4.2 AVRCOMB, NEW10.
p$data<-df_with_tests_new|>filter(attribute%in%"AVRCOMB")
p
p$data<-df_with_tests_new|>filter(attribute%in%"NEW10")
p
df_with_tests|>dplyr::filter(attribute%in%colNew[-length(colNew)])|>dplyr::select(-PPIND, -value)|>slice_head(n=2)|>rmarkdown::paged_table()2.4.2.1 Наблюдения:
Можем сказать, что public и private для AVRCOMB и NEW10 имеют разные матожидания.
В случае гипотезы о равенстве дисперсий, критерий Левена (Levene) и критерий Фишера (var) имеют разные мощности. Критерий Фишера используется, если сравниваемые выборки имеют нормальное распределение. Поскольку группа public признака NEW10 отвергла гипотезу о нормальном распределении, применять критерий Фишера нельзя.
Отвергаем равенство дисперсий для AVRCOMB, NEW10.
2.4.3 BOOK, log-ADD_FEE.
p$data<-df_with_tests_fin|>filter(attribute%in%"BOOK")
p
p$data<-df_with_tests_fin|>filter(attribute%in%"log-ADD_FEE")
p
df_with_tests|>dplyr::filter(attribute%in%c("log-ADD_FEE", "BOOK"))|>dplyr::select(-PPIND, -value)|>head(length(colFin)-1)|>rmarkdown::paged_table()2.4.3.1 Наблюдения:
Исходя из p_value не отвергаем, что public и private имеют равные матожидания. Тест Вилкоксона также характеризует их однородность.
2.4.3.2 Задача 1.1 из .tsk
Проверим log-ADD_FEE и BOOK друг с другом
compare_column_tests <- function(column_1,
column_2,
alpha = 0.05) {
levene_test<-leveneTest(column_2, column_1)
var_test <- var.test(column_1, column_2)
wilcox_test <- wilcox.test(column_1, column_2)
t_test_equal <- t.test(column_1, column_2, var.equal = TRUE)
t_test_unequal <- t.test(column_1, column_2, var.equal = FALSE)
results <- data.frame(
t_test_unequal_p_value = t_test_unequal$p.value,
t_test_equal_p_value = t_test_equal$p.value,
levene_test_p_value = levene_test$`Pr(>F)`[1],
var_test_p_value = var_test$p.value,
wilcox_test_p_value = wilcox_test$p.value
)
}af_pr<-col_I_sn$ADD_FEE[col_I_sn$PPIND=="private"]
af_pu<-col_I_sn$ADD_FEE[col_I_sn$PPIND=="public"]
b_pr<-col_split_mixed$private$BOOK
b_pu<-col_split_mixed$public$BOOK
print(compare_column_tests(af_pr, b_pr))## t_test_unequal_p_value t_test_equal_p_value levene_test_p_value
## 1 0.005903872 0.003538701 2.925758e-75
## var_test_p_value wilcox_test_p_value
## 1 1.42501e-07 1.678377e-05print(compare_column_tests(af_pu, b_pu))## t_test_unequal_p_value t_test_equal_p_value levene_test_p_value
## 1 0.1706573 0.1115642 0.9495775
## var_test_p_value wilcox_test_p_value
## 1 0 0.00060916512.4.3.2.1 Наблюдения:
В общем случае тесты показывают, что разница есть в случае общественных вузов, но для частных вузов не отвергаем равенство дисперсий (из levene’s test) и матожиданий.
2.4.4 GRADUAT, PH_D, SF_RATIO
p$data<-df_with_tests_stud
p
df_with_tests|>dplyr::filter(attribute%in%colStud[-length(colStud)])|>dplyr::select(-PPIND, -value)|>head(length(colStud)-1)|>rmarkdown::paged_table()2.4.4.1 Наблюдения:
Отвергаем равенство матожиданий для GRADUAT, SF_RATIO, PH_D. Не отвергаем равенство дисперсий GRADUAT, SF_RATIO, PH_D. — ### OUT_STAT, R_B_COST
p$data<-df_with_tests_cost
p
df_with_tests|>dplyr::filter(attribute%in%colCost[-length(colCost)])|>dplyr::select(-PPIND, -value)|>head(length(colCost)-2)|>rmarkdown::paged_table()2.4.4.2 Наблюдения:
Отвергаем равенство матожиданий для OUT_STAT, R_B_COST. Не отвергаем равенство дисперсий для R_B_COST.
2.5 Далее, объясняете, в каких случаях (распределение далеко от нормального, могут быть выделяющиеся наблюдения, нас интересует другая характеристика положения) t-критерий не удовлетворителен и нужно переходить к непараметрическим критериям. Рассказываете, какой из непараметрических критериев является аналогом t-критерия, как он строится и против какой альтернативы мощный. Вы уже догадались, что это критерий Манна-Уитни, он же критерий Вилкоксона.
t-критерий неудовлетворителен в нескольких ситуациях:
- Малый объём выборки. Если форма распределения признака неизвестна или не нормальна, то t-критерий асимптотически точный.
- Наличие выбросов. t-критерий неусточив по отношению к выбросам.
- Другие цели эксперимента. t-критерий мощный против альтернатив о неравенстве матожиданий, его нельзя применить для построения других гипотез, например, гипотезы о равенстве дисперсий.
В качестве замены t-критерия можно воспользоваться непараметрическим критерием Манна-Уитни (Вилкоксона).
2.6 Смотрите на результаты применения критерия Манна-Уитни, сравниваете с результатами применения t-критерия. Проводите сравнительный анализ критериев с теоретической точки зрения (чем один лучше другого и чем хуже).
df_tests|>dplyr::select(attribute,t_test_unequal_p_value, t_test_equal_p_value, wilcox_test_p_value)|>rmarkdown::paged_table()2.6.1 Наблюдения:
В случаях, где Манн-Уитни мощнее t-test: BOOK, R_B_COST, OUT_STAT, SF_RATIO, GRADUAT; причиной этого могут быть выбросы. В случаях, где Манн-Уитни не мощнее t-test: log-ADD_FEE, NEW10, AVRCOMB, PH_D; могут отсутствовать выбросы, из-за чего t-test оказывается мощнее.
2.6.2 Теоретическое сравнение:
Как непараметрический критерий, тест Вилкоксона уступает t-тесту в мощности против альтернативы о неравенстве матожиданий. С другой стороны, тест Вилкоксона устойчив к выбросам, мощный против альтернативы о неравенстве распределения двух случайных величин. Если известно, что две случайные величины имеют одинаковые формы распределения, то тест Вилкоксона автоматически ставит гипотезу о равенстве матожиданий.
2.7 Далее, переходите к критериям, которые умеют сравнивать не только характеристики положения, но и формы распределений (это критерий Колмогорова-Смирнова, например). Для каждого критерия (включая критерий Манна-Уитни), нужно уметь объяснять, что означают столбцы в таблицах результатов критериев. Также, при разных результатах проверки гипотезы о равенстве распределений нужно объяснять, почему один критерий, например, не отверг гипотезу, а другой – отверг.
2.7.1 Сравнение форм распределения public и private
df_tests|>dplyr::select(attribute, wilcox_test_p_value, ks_test_p_value)|>rmarkdown::paged_table()2.7.1.1 Наблюдения:
Столбец wilcox_test_p_value - p-value из критерия Вилкоксона. Столбец ks_test_p_value - p-value из критерия Колмогорова-Смирнова. Критерии Вилкоксона и Колмогорова-Смирнова мощны против альтернативы о неравенстве распределения. Критерии не отвергают равенство форм распределений для log-ADD_FEE и BOOK, но отвергают для остальных. Среди рассматриваемых признаков нет случаев, когда один критерий отвергает гипотезу, а другой - не отвергает.
3 Об анализе зависимостей
3.1 Вспомните, какие бывают виды зависимостей и чем они измеряются, по каким формулам. Посмотрите на основе pairs plot, какие зависимости у вас в данных. Не забудьте, что при неоднородных данных изучать зависимости имеет смысл только внутри групп по-отдельности.
col_split_mixed<-split(col_I_mixed, col_I_mixed$PPIND)
names(col_split_mixed)<-c("private", "public")
col_split_mixed<-lapply(col_split_mixed, function(x){
dplyr::select(x, -PPIND)
})
colMixed<-names(col_I_mixed)
colMixed<-colMixed[-length(colMixed)]
col_split_mixed$private|>rmarkdown::paged_table()col_split_mixed$public|>rmarkdown::paged_table()3.1.1 Private:
ggpairs(col_split_mixed$private, aes(alpha = 0.5),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))
3.1.1.1 Наблюдения:
Можно заметить, что NEW10/AVRCOMB разбивается на две части (поскольку AVRCOMB имеет много пропусков и сам AVRCOMB подразумевает разделение на подгруппы).
col_I_sn[col_I_sn$PPIND=="private",]|>filter(NEW10<=50 & !is.na(AVRCOMB))|>rmarkdown::paged_table()col_I_sn[col_I_sn$PPIND=="private",]|>filter(NEW10>50 & !is.na(AVRCOMB))|>rmarkdown::paged_table()Среди частных университетов только один “Delaware University” имеет расхождения между IN_STATE и OUT_STAT. Сейчас идёт такое ценообразование.
3.1.1.1.1 University of Delaware.
University of Delaware’s tuition is $15,410 for in-state and $37,930 for out-of-state.
col_I_sn[col_I_sn$IN_STATE!=col_I_sn$OUT_STAT & col_I_sn$PPIND=="private",]3.1.2 Public:
ggpairs(col_split_mixed$public, aes(alpha = 0.5),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))
3.1.2.1 Наблюдения
На pairs plot можно заметить выбросы в признаках log-ADD_FEE, OUT_STAT, NEW10, которые влияют на коэффициент корреляции Пирсона.
3.2 Начинать нужно с анализа линейных зависимостей. На основе коэффициента корреляции Пирсона нужно проинтерпретировать значимые зависимости. При наличие в данных пропусков обратите внимание на выбор между casewise and pairwise MD deletion (в чем разница, какие недостатки и достоинства у этих вариантов?).
Учитывая ограниченный объём данных, тем более после разделения по группам (например, наименьшая в private размера 29 без NA), следует воспользоваться pairwise. Более того, из-за признака AVRCOMB использование casewise неявно рассматривает корреляцию в другой подгруппе. Минусом pairwise является то, что из-за попарной проверки могут учитываться данные, пропущенные в других парах.
3.3 Затем можно переходить к ранговым коэффициентам корреляции. Расскажите, при каких условиях коэффициенты корреляции Пирсона и Спирмена примерно равны. Приведите примеры, когда один из них больше другого и наоборот. Сравните результаты на ваших данных. Если при сравнении будут найдены заметные различия в результатах, то попробуйте объяснить причину.
compute_corr <- function(group_data) {
corr_pairwise_pearson <- cor(group_data, use = "pairwise.complete.obs", method = "pearson")
corr_pairwise_spearman <- cor(group_data, use = "pairwise.complete.obs", method = "spearman")
return(list(pairwise_pearson = corr_pairwise_pearson, pairwise_spearman = corr_pairwise_spearman))
}
make_corr_plots <- function(corr_list, group) {
p1 <- ggcorrplot(corr_list$pairwise_pearson,
lab = TRUE,
lab_size = 1.9,
colors = c("blue", "white", "red"),
title = paste("Pairwise Cor Pearson Matrix for Group", group),
ggtheme = theme_minimal() +
theme(plot.title = element_text(size = 6)) # Reduced title size
)
p2 <- ggcorrplot(corr_list$pairwise_spearman,
lab = TRUE,
lab_size = 1.9,
colors = c("blue", "white", "red"),
title = paste("Pairwise Cor Spearman Matrix for Group", group),
ggtheme = theme_minimal() +
theme(plot.title = element_text(size = 6)) # Reduced title size
)
grid.arrange(p1, p2, ncol = 2, widths = c(4, 4), # Adjust widths here
padding = unit(1, "cm")) # Add padding for extra space
}library(ggcorrplot)
for (group in names(col_split_mixed)) {
group_data <- col_split_mixed[[group]]
cor_matrix <- compute_corr(group_data)
plot_pearson_spearman <- make_corr_plots(cor_matrix, group)
print(plot_pearson_spearman)
}
## TableGrob (1 x 2) "arrange": 2 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## TableGrob (1 x 2) "arrange": 2 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]3.3.1 Наблюдения, различия между коэффициентом корреляции Пирсона и Спирмена:
3.3.1.1 Private:
Пирсон мощнее Спирмена:
- OUT_STAT/R_B_COST - почти линейный график и выбросы,
- SF_RATIO/BOOK - выбросы при данных, похожих на облако.
- R_B_COST/PH_D - выбросы усиливают линейную зависимость,
- R_B_COST/GRADUAT - выбросы усиливают линейную зависимость.
Спирмен мощнее Пирсона:
- OUT_STAT/NEW10 - нелинейные данные,
- R_B_COST/BOOK - выбросы,
- log-ADD_FEE/OUT_STAT - выбросы,
- BOOK/GRADUAT - выбросы.
3.3.1.2 Public:
Пирсон мощнее Спирмена:
- R_B_COST/log-ADD_FEE - выбросы,
- NEW10/log-ADD_FEE - выбросы усиливают линейную зависимость.
Спирмен мощнее Пирсона:
- OUT_STAT/log-ADD_FEE - выбросы,
- AVRCOMB/NEW10 - выбросы уменьшают коэффициент корреляции Пирсона.
3.3.2 Теоретическое объяснение (сравнение коэффициента корреляции Пирсона и Спирмена):
В тех ячейках, где совпадают или близки коэффициенты корреляции Пирсона и Спирмена почти нет выбросов и у scatterplot эллиптическая структура, если Пирсон значительно больше Спирмена - есть заметные выбросы. Спирмен может быть больше Пирсона, если выбросы имеются, но их положение не очень значимо относительно других наблюдений.
Значительные корреляции:
3.3.3 Private:
- Выраженный блок GRADUAT, PH_D, SF_RATIO, AVRCOMB, NEW10, OUT_STAT.
3.3.4 Public:
- Блок NEW10, AVRCOMB, PH_D, GRADUAT;
- Пара NEW10, R_B_COST - нельзя сказать (по смыслу) о зависимости переменных. Могут зависеть от внешней переменной (скрытого фактора).
3.4 Проинтерпретируйте найденные корреляции - можно ли сказать, что является причиной, что следствием. Если есть какая-то другая причина, которая влияет одновременно на оба признака (скрытый фактор), то попробуйте убрать его влияние с помощью частных корреляций.
3.4.1 Private:
3.4.2 Public:
4 Регрессионный анализ
Задача: провести регрессию качества поступающих студентов (NEW10) на остальные признаки. С помощью пошаговой регрессии уменьшить число признаков
library(lm.beta)4.1 Подготовка данных
4.1.1 Выбор признаков
Рассматриваемые наблюдения: 1. PPIND (для группировки) 2. ADD_FEE - дополнительные выплаты 3. BOOK - плата за бронирование 4. NEW10 (зависимая переменная, также NEW25) - студентов из топ 10% своей старшей школы 5. AVRCOMB - средний комбинированный балл 6. SF_RATIO - Student/Faculty ratio 7. PH_D - Pct. of faculty with Ph.D.’s 8. GRADUAT - Graduation rate 9. R_B_COST - Room and board costs 10. OUT_STAT - Out-of-state tuition
colNames<-c("...1", "ADD_FEE", "BOOK", "R_B_COST","OUT_STAT","NEW10","AVRCOMB","SF_RATIO","PH_D","GRADUAT","PPIND")
col_I_regr<-col_I_sn[c(colNames)]
col_I_regr<-col_I_regr|>mutate(ADD_FEE = log(ADD_FEE))|>rename(log_ADD_FEE=ADD_FEE)
col_I_regr_split<-split(col_I_regr, col_I_regr$PPIND)# Автор: -Николай-
library(GGally)
library(ggplot2)
library(dplyr)
library(rlang)##
## Присоединяю пакет: 'rlang'## Следующие объекты скрыты от 'package:purrr':
##
## %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
## flatten_raw, invoke, spliceplot_numeric_pairs <- function(df, title = NULL, highlight = NULL) {
# Выбираем имена всех numeric-столбцов
numeric_cols <- df %>% dplyr::select(dplyr::where(is.numeric)) %>% names()
# Функция для нижних панелей с возможностью выделения точек разными цветами
custom_points <- function(data, mapping, ...) {
# Отрисовываем базовые точки синим цветом
p <- ggplot(data, mapping) +
geom_point(alpha = 0.6, size = 2, color = "blue", ...)
if (!is.null(highlight)) {
# Если highlight не является списком, превращаем его в список с одним элементом и цветом "red"
if (!is.list(highlight) || inherits(highlight, "data.frame")) {
highlight <- list(red = highlight)
}
# Получаем имена переменных из mapping с помощью as_label()
xvar <- as_label(mapping$x)
yvar <- as_label(mapping$y)
# Для каждого набора точек в списке highlight
for (col in names(highlight)) {
hl_df <- highlight[[col]]
# Проверяем, что hl_df содержит нужные переменные
if(all(c(xvar, yvar) %in% names(hl_df))) {
# Находим совпадающие строки между данными панели и hl_df
data_highlight <- semi_join(data, hl_df, by = c(xvar, yvar))
# Накладываем выделенные точки нужного цвета
p <- p + geom_point(data = data_highlight, mapping = mapping,
alpha = 0.6, size = 2, color = col, ...)
}
}
}
p
}
# Строим матрицу парных графиков
ggpairs(df,
columns = numeric_cols,
columnLabels = numeric_cols,
title = title,
upper = list(continuous = wrap("cor", size = 3)),
lower = list(continuous = custom_points),
diag = list(continuous = wrap("barDiag", bins = 10, fill = "lightblue", color = "black"))
) +
theme_bw() +
theme(
strip.text.x = element_text(size = 7, angle = 0, hjust = 0),
strip.text.y = element_text(size = 8, angle = 0, hjust = 0),
axis.text = element_text(size = 6),
plot.margin = unit(c(1, 2, 1, 1), "cm")
)
}4.2 Анализ частных вузов
pr_data<-col_I_regr_split$private
print(paste0("Число наблюдений: ", dim(pr_data)[1]))## [1] "Число наблюдений: 53"colSums(is.na(pr_data))## ...1 log_ADD_FEE BOOK R_B_COST OUT_STAT NEW10
## 0 7 0 0 0 1
## AVRCOMB SF_RATIO PH_D GRADUAT PPIND
## 24 0 4 1 0Много пропусков в AVRCOMB
plot_numeric_pairs(pr_data, "Частные вузы")
4.2.1 Выбросы “на глаз”
outlier_pr_book_eye<-pr_data[pr_data$BOOK>875,]
outlier_pr_out_eye<-pr_data[pr_data$OUT_STAT<5000,]
outlier_pr_book_eyeoutlier_pr_out_eyeplot_numeric_pairs(pr_data, "Частные вузы", highlight = list(red=outlier_pr_book_eye, black=outlier_pr_out_eye))
4.2.2 Построение линейной регрессии
Здесь условия из .tsk файла
Стоит помнить, что у AVRCOMB много NA. В такой модели не будет выбросов “на глаз”
model_pr_new <- lm(NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT,
data = pr_data)|>lm.beta()4.2.3 Интерпретация результатов
summary(model_pr_new)##
## Call:
## lm(formula = NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST +
## PH_D + SF_RATIO + GRADUAT + OUT_STAT, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.4148 -3.8409 0.0271 2.9131 12.4176
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -78.090327 NA 51.790742 -1.508 0.1575
## AVRCOMB 0.117625 0.544789 0.042417 2.773 0.0169 *
## BOOK -0.006699 -0.029617 0.014003 -0.478 0.6410
## log_ADD_FEE 3.200735 0.098598 2.219759 1.442 0.1749
## R_B_COST -0.003437 -0.117762 0.002499 -1.376 0.1941
## PH_D -0.492675 -0.126378 0.496416 -0.992 0.3406
## SF_RATIO -1.353492 -0.205756 0.594526 -2.277 0.0419 *
## GRADUAT 0.306142 0.178421 0.258574 1.184 0.2594
## OUT_STAT 0.002354 0.257544 0.001391 1.693 0.1163
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.797 on 12 degrees of freedom
## (32 пропущенных наблюдений удалены)
## Multiple R-squared: 0.9646, Adjusted R-squared: 0.9411
## F-statistic: 40.92 on 8 and 12 DF, p-value: 1.496e-074.2.4 Построение доверительных эллипсоидов
# Автор: -Вадим-
create_ellipse <- function(center, cov_matrix, level = 0.95, n = 100) {
# Уровень доверия преобразуем в квантиль хи-квадрат
chi_sq <- qchisq(level, df = 2)
# Вычисляем собственные значения и собственные векторы
eigen_res <- eigen(cov_matrix)
# Углы для точек эллипса
angles <- seq(0, 2 * pi, length.out = n)
# Радиусы эллипса (основаны на собственных значениях)
radii <- sqrt(chi_sq * eigen_res$values)
# Создаем эллипс в стандартной ориентации
ellipse_points <- cbind(radii[1] * cos(angles), radii[2] * sin(angles))
# Поворачиваем эллипс в соответствии с собственными векторами
ellipse_points <- t(eigen_res$vectors %*% t(ellipse_points))
# Добавляем смещение
ellipse_points <- sweep(ellipse_points, 2, center, "+")
ellipse_points <- as.data.frame(ellipse_points)
colnames(ellipse_points) <- c("x", "y") # Assign column names
return(ellipse_points)
}
plot_confidence_ellipse <- function(model, which.coeff = c(2, 4), level = 0.95,
xlab = NULL, ylab = NULL, main = "Доверительный эллипс",
margin = 0.1) {
# Проверка на то, что переданы именно 2 коэффициента
if (length(which.coeff) != 2) {
stop("which.coeff должен содержать ровно 2 индекса коэффициентов")
}
# Получение необходимых данных из модели
model_summary <- summary(model)
coef_estimates <- coef(model)
сov_matrix <- model_summary$cov.unscaled
# Проверка допустимости индексов коэффициентов
if (any(which.coeff > length(coef_estimates))) {
stop("Индекс коэффициента превышает их количество в модели")
}
# Создание эллипса
ellipse_points <- create_ellipse(
center = coef_estimates[which.coeff],
cov_matrix = сov_matrix[which.coeff, which.coeff],
level = level
)
# Определение границ графика с учетом отступов
x_min <- min(ellipse_points$x)
x_max <- max(ellipse_points$x)
y_min <- min(ellipse_points$y)
y_max <- max(ellipse_points$y)
x_range <- c(x_min - margin * (x_max - x_min), x_max + margin * (x_max - x_min))
y_range <- c(y_min - margin * (y_max - y_min), y_max + margin * (y_max - y_min))
# Если не указаны названия осей, используем имена коэффициентов
if (is.null(xlab)) {
xlab <- names(coef_estimates)[which.coeff[1]]
}
if (is.null(ylab)) {
ylab <- names(coef_estimates)[which.coeff[2]]
}
# Создаем data.frame с центральной точкой
center_point <- data.frame(
x = coef_estimates[which.coeff[1]],
y = coef_estimates[which.coeff[2]]
)
# Создаем объект ggplot
p <- ggplot2::ggplot() +
# Добавляем эллипс
ggplot2::geom_path(data = ellipse_points, ggplot2::aes(x = x, y = y), color = "blue") +
# Добавляем точку с оценкой коэффициента
ggplot2::geom_point(data = center_point, ggplot2::aes(x = x, y = y), color = "red", size = 3) +
# Добавляем вертикальную и горизонтальную линии через центр
ggplot2::geom_vline(xintercept = center_point$x, linetype = "dashed", color = "gray") +
ggplot2::geom_hline(yintercept = center_point$y, linetype = "dashed", color = "gray") +
# Настраиваем оси и заголовок
ggplot2::xlim(x_range) +
ggplot2::ylim(y_range) +
ggplot2::labs(
x = xlab,
y = ylab,
title = main
) +
ggplot2::theme_light()
return(p)
}plot_confidence_ellipse(model_pr_new, c(3, 2), level=0.95)
На примере этой пары: AVRCOMB и SF_RATIO значимы
4.2.5 Влияние корреляции на качество оценки
На примере, PH_D и OUT_STAT
phcor<-cor(pr_data$PH_D, pr_data$OUT_STAT, use = "complete.obs")
y <- pr_data$NEW10
X <- pr_data[, c("PH_D", "OUT_STAT")]
model_cor <- lm(NEW10 ~ PH_D + OUT_STAT, data = pr_data)
sum_cor<-summary(model_cor)
sum_cor##
## Call:
## lm(formula = NEW10 ~ PH_D + OUT_STAT, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.696 -13.720 0.979 15.266 34.815
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.009e+02 3.392e+01 -2.975 0.00470 **
## PH_D 1.551e+00 4.763e-01 3.257 0.00214 **
## OUT_STAT 1.198e-03 9.551e-04 1.255 0.21608
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.96 on 45 degrees of freedom
## (5 пропущенных наблюдений удалены)
## Multiple R-squared: 0.4428, Adjusted R-squared: 0.4181
## F-statistic: 17.88 on 2 and 45 DF, p-value: 1.926e-06# Оцениваем дисперсию ошибок
sigma_s_squared <- (sum_cor$sigma/sd(y, na.rm=TRUE))^2 # оценка дисперсии ошибок
# Размер выборки
n <- nobs(model_cor)
R_XX<-cor(X, use="complete.obs")
print("Inverse R_XX")## [1] "Inverse R_XX"solve(R_XX)## PH_D OUT_STAT
## PH_D 2.088331 -1.507579
## OUT_STAT -1.507579 2.088331cov_beta<-sigma_s_squared/n*solve(R_XX)
print("Cov_beta:")## [1] "Cov_beta:"cov_beta## PH_D OUT_STAT
## PH_D 0.02360939 -0.01704377
## OUT_STAT -0.01704377 0.02360939print("Determinant")## [1] "Determinant"detcoef<-1/(1-phcor^2)
detcoef## [1] 2.088331print("-rho/(1-rho^2)")## [1] "-rho/(1-rho^2)"-phcor*detcoef## [1] -1.507579phcor## [1] 0.72190634.2.6 Множественная корреляция
# Автор: -Евгений-
# Вычисление множественной корреляции
calc_multicollinearity <- function(var, data) {
X <- as.matrix(data[, setdiff(names(data), var)]) # Матрица остальных регрессоров
y <- as.matrix(data[[var]]) # Текущий регрессор как зависимая переменная
model <- lm(y ~ X)
sqrt(summary(model)$r.squared)
}
calc_partial_correlation <- function(var, data) {
# Остатки регрессии текущего регрессора на остальные регрессоры
data<-data|>na.omit()
model_X <- lm(data[[var]] ~ ., data = data[, !names(data) %in% c(var, "NEW10")])
res_X <- residuals(model_X)
# Остатки регрессии Y на остальные регрессоры
model_Y <- lm(NEW10 ~ ., data = data[, !names(data) %in% var])
res_Y <- residuals(model_Y)
cor(res_X, res_Y, use = "complete.obs") # Корреляция остатков
}
regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)
print(redundancy_table)## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.6803712 0.3842868
## BOOK BOOK 0.4955889 -0.1368050
## R_B_COST R_B_COST 0.8079355 -0.3690491
## OUT_STAT OUT_STAT 0.9472357 0.4390291
## AVRCOMB AVRCOMB 0.9764603 0.6249356
## SF_RATIO SF_RATIO 0.8649231 -0.5492085
## PH_D PH_D 0.9121799 -0.2754193
## GRADUAT GRADUAT 0.9401219 0.3234121Сильно положительные множественные корреляции вместе с частными корреляциями, следует сравнить по отдельности
regr_high <- c("R_B_COST", "OUT_STAT", "AVRCOMB", "SF_RATIO", "PH_D", "GRADUAT")
regr_fullhigh<-c(regr_high, "NEW10")
redund_high <- data.frame(
Регрессор = regr_high,
Множественная_корреляция = sapply(regr_high, calc_multicollinearity, pr_data[, regr_fullhigh]),
Частная_корреляция = sapply(regr_high, calc_partial_correlation, pr_data[,regr_fullhigh])
)
print(redund_high)## Регрессор Множественная_корреляция Частная_корреляция
## R_B_COST R_B_COST 0.6621454 -0.1634685
## OUT_STAT OUT_STAT 0.9214250 0.3142483
## AVRCOMB AVRCOMB 0.9687653 0.5706967
## SF_RATIO SF_RATIO 0.8369497 -0.3982266
## PH_D PH_D 0.8830917 -0.1158525
## GRADUAT GRADUAT 0.9332980 0.3403546Видно, что PH_D имеет высокую множественную корреляцию и низкую частную корреляцию, поэтому выбросим её
model_pr_noPhD<-lm(NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + SF_RATIO + GRADUAT + OUT_STAT, data=pr_data)|>lm.beta()
summary(model_pr_noPhD)##
## Call:
## lm(formula = NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST +
## SF_RATIO + GRADUAT + OUT_STAT, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.9908 -2.5824 0.6409 3.1699 11.8038
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -1.156e+02 NA 2.833e+01 -4.080 0.000872 ***
## AVRCOMB 1.291e-01 6.067e-01 2.746e-02 4.703 0.000240 ***
## BOOK -4.336e-03 -1.885e-02 1.215e-02 -0.357 0.725883
## log_ADD_FEE 2.568e+00 7.431e-02 2.010e+00 1.277 0.219712
## R_B_COST -3.258e-03 -1.042e-01 2.189e-03 -1.488 0.156128
## SF_RATIO -1.300e+00 -1.876e-01 5.055e-01 -2.572 0.020479 *
## GRADUAT 1.590e-01 8.799e-02 1.657e-01 0.960 0.351347
## OUT_STAT 1.837e-03 1.893e-01 1.079e-03 1.704 0.107803
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.753 on 16 degrees of freedom
## (29 пропущенных наблюдений удалены)
## Multiple R-squared: 0.9616, Adjusted R-squared: 0.9448
## F-statistic: 57.28 on 7 and 16 DF, p-value: 3.83e-10AIC(model_pr_noPhD)## [1] 168.0601regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10", "PH_D"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)
print(redundancy_table)## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.5973852 0.30419786
## BOOK BOOK 0.3840774 -0.08885404
## R_B_COST R_B_COST 0.7551732 -0.34870910
## OUT_STAT OUT_STAT 0.9141040 0.39183509
## AVRCOMB AVRCOMB 0.9692789 0.76171345
## SF_RATIO SF_RATIO 0.8252002 -0.54081295
## GRADUAT GRADUAT 0.8544490 0.23336928Попробуем убрать BOOK (из-за низкой частной корреляции) и GRADUAT
model_pr1<-lm(NEW10 ~ AVRCOMB + log_ADD_FEE + R_B_COST + SF_RATIO + OUT_STAT, data=pr_data)|>lm.beta()
summary(model_pr1)##
## Call:
## lm(formula = NEW10 ~ AVRCOMB + log_ADD_FEE + R_B_COST + SF_RATIO +
## OUT_STAT, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.4607 -1.9482 -0.7267 2.3280 13.4017
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -1.185e+02 NA 2.745e+01 -4.315 0.000417 ***
## AVRCOMB 1.372e-01 6.446e-01 2.323e-02 5.907 1.36e-05 ***
## log_ADD_FEE 2.961e+00 8.567e-02 1.880e+00 1.575 0.132718
## R_B_COST -3.886e-03 -1.243e-01 2.040e-03 -1.905 0.072923 .
## SF_RATIO -1.336e+00 -1.927e-01 4.914e-01 -2.718 0.014095 *
## OUT_STAT 2.164e-03 2.229e-01 9.870e-04 2.192 0.041756 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.583 on 18 degrees of freedom
## (29 пропущенных наблюдений удалены)
## Multiple R-squared: 0.959, Adjusted R-squared: 0.9476
## F-statistic: 84.15 on 5 and 18 DF, p-value: 7.676e-12AIC(model_pr1)## [1] 165.6631regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10", "PH_D", "BOOK", "GRADUAT"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)
print(redundancy_table)## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.5684235 0.3479827
## R_B_COST R_B_COST 0.7447506 -0.4095647
## OUT_STAT OUT_STAT 0.9088596 0.4590401
## AVRCOMB AVRCOMB 0.9668878 0.8122173
## SF_RATIO SF_RATIO 0.8236772 -0.5394659AIC(model_pr_new)## [1] 148.3376AIC(model_pr_noPhD)## [1] 168.0601AIC(model_pr1)## [1] 165.66314.2.7 OLS. Пошаговый AIC.
library(olsrr)##
## Присоединяю пакет: 'olsrr'## Следующий объект скрыт от 'package:MASS':
##
## cement## Следующий объект скрыт от 'package:datasets':
##
## riversbackward_aic<-ols_step_backward_aic(model_pr_new, progress=TRUE, details=TRUE)## Backward Elimination Method
## ---------------------------
##
## Candidate Terms:
##
## 1. AVRCOMB
## 2. BOOK
## 3. log_ADD_FEE
## 4. R_B_COST
## 5. PH_D
## 6. SF_RATIO
## 7. GRADUAT
## 8. OUT_STAT
##
##
## Step => 0
## Model => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 148.3376
##
## Initiating stepwise selection...
##
## Table: Removing Existing Variables
## -----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -----------------------------------------------------------------------
## BOOK 1 146.734 156.135 97.586 0.96397 0.94457
## PH_D 1 147.994 157.395 97.513 0.96174 0.94114
## GRADUAT 1 148.658 158.058 97.491 0.96051 0.93925
## R_B_COST 1 149.412 158.813 97.480 0.95907 0.93703
## log_ADD_FEE 1 149.693 159.094 97.480 0.95852 0.93618
## OUT_STAT 1 150.834 160.235 97.501 0.95620 0.93262
## SF_RATIO 1 153.877 163.277 97.737 0.94937 0.92211
## AVRCOMB 1 156.737 166.137 98.197 0.94199 0.91075
## -----------------------------------------------------------------------
##
## Step => 1
## Removed => BOOK
## Model => NEW10 ~ AVRCOMB + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 146.7343
##
## Table: Removing Existing Variables
## -----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -----------------------------------------------------------------------
## PH_D 1 147.170 155.526 94.423 0.95954 0.94220
## GRADUAT 1 147.861 156.217 94.521 0.95818 0.94026
## log_ADD_FEE 1 147.962 156.319 94.536 0.95798 0.93997
## R_B_COST 1 148.265 156.622 94.583 0.95737 0.93910
## OUT_STAT 1 150.933 159.289 95.079 0.95160 0.93085
## SF_RATIO 1 153.017 161.373 95.578 0.94655 0.92364
## AVRCOMB 1 155.056 163.412 96.162 0.94110 0.91585
## -----------------------------------------------------------------------
##
##
## No more variables to be removed.
##
## Variables Removed:
##
## => BOOKforward_aic<-ols_step_forward_aic(model_pr_new, progress=TRUE, details=TRUE)## Forward Selection Method
## ------------------------
##
## Candidate Terms:
##
## 1. AVRCOMB
## 2. BOOK
## 3. log_ADD_FEE
## 4. R_B_COST
## 5. PH_D
## 6. SF_RATIO
## 7. GRADUAT
## 8. OUT_STAT
##
##
## Step => 0
## Model => NEW10 ~ 1
## AIC => 202.5251
##
## Initiating stepwise selection...
##
## Table: Adding New Variables
## -------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## AVRCOMB 1 147.617 150.751 87.731 0.93346 0.92995
## OUT_STAT 1 171.142 174.276 107.789 0.79601 0.78527
## GRADUAT 1 176.803 179.937 112.954 0.73289 0.71883
## PH_D 1 186.429 189.563 121.963 0.57756 0.55533
## SF_RATIO 1 189.574 192.707 124.955 0.50933 0.48350
## log_ADD_FEE 1 202.476 205.609 137.417 0.09297 0.04523
## R_B_COST 1 203.410 206.543 138.328 0.05173 0.00182
## BOOK 1 204.456 207.590 139.350 0.00327 -0.04919
## -------------------------------------------------------------------------
##
## Step => 1
## Added => AVRCOMB
## Model => NEW10 ~ AVRCOMB
## AIC => 147.617
##
## Table: Adding New Variables
## -----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -----------------------------------------------------------------------
## SF_RATIO 1 145.049 149.227 86.340 0.94647 0.94052
## OUT_STAT 1 147.493 151.671 88.069 0.93986 0.93318
## BOOK 1 147.975 152.153 88.412 0.93846 0.93162
## log_ADD_FEE 1 148.581 152.759 88.846 0.93666 0.92962
## PH_D 1 149.095 153.274 89.215 0.93509 0.92788
## GRADUAT 1 149.166 153.344 89.265 0.93487 0.92764
## R_B_COST 1 149.275 153.453 89.344 0.93453 0.92726
## -----------------------------------------------------------------------
##
## Step => 2
## Added => SF_RATIO
## Model => NEW10 ~ AVRCOMB + SF_RATIO
## AIC => 145.0489
##
## Table: Adding New Variables
## -----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -----------------------------------------------------------------------
## BOOK 1 145.266 150.488 87.603 0.95082 0.94215
## OUT_STAT 1 145.397 150.620 87.679 0.95052 0.94178
## GRADUAT 1 146.151 151.373 88.121 0.94871 0.93966
## log_ADD_FEE 1 146.479 151.701 88.314 0.94790 0.93871
## R_B_COST 1 146.891 152.113 88.557 0.94687 0.93749
## PH_D 1 146.968 152.190 88.603 0.94667 0.93726
## -----------------------------------------------------------------------
##
##
## No more variables to be added.
##
## Variables Selected:
##
## => AVRCOMB
## => SF_RATIOПолучили, что forward_aic имеет лучшее значение, чем backward_aic, то есть, является более хорошей линейной моделью. К тому же, у forward_aic меньше признаков.
4.2.7.1 Forward
forward_aic##
##
## Stepwise Summary
## -------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## 0 Base Model 202.525 204.614 139.293 0.00000 0.00000
## 1 AVRCOMB 147.617 150.751 87.731 0.93346 0.92995
## 2 SF_RATIO 145.049 149.227 86.340 0.94647 0.94052
## -------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ---------------------------------------------------------------
## R 0.973 RMSE 6.323
## R-Squared 0.946 MSE 39.975
## Adj. R-Squared 0.941 Coef. Var 11.813
## Pred R-Squared 0.925 AIC 145.049
## MAE 4.914 SBC 149.227
## ---------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## ----------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ----------------------------------------------------------------------
## Regression 14841.759 2 7420.879 159.117 0.0000
## Residual 839.479 18 46.638
## Total 15681.238 20
## ----------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------------
## (Intercept) -150.963 21.538 -7.009 0.000 -196.212 -105.714
## AVRCOMB 0.188 0.015 0.869 12.124 0.000 0.155 0.220
## SF_RATIO -0.986 0.471 -0.150 -2.091 0.051 -1.977 0.004
## ----------------------------------------------------------------------------------------------4.2.7.2 Backward
backward_aic##
##
## Stepwise Summary
## -------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## 0 Full Model 148.338 158.783 101.117 0.96464 0.94107
## 1 BOOK 146.734 156.135 96.228 0.96397 0.94457
## -------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ---------------------------------------------------------------
## R 0.982 RMSE 5.187
## R-Squared 0.964 MSE 26.905
## Adj. R-Squared 0.945 Coef. Var 11.404
## Pred R-Squared 0.901 AIC 146.734
## MAE 4.081 SBC 156.135
## ---------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## ---------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ---------------------------------------------------------------------
## Regression 15116.223 7 2159.460 49.685 0.0000
## Residual 565.015 13 43.463
## Total 15681.238 20
## ---------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------------
## (Intercept) -71.623 48.490 -1.477 0.163 -176.379 33.133
## AVRCOMB 0.110 0.038 0.511 2.873 0.013 0.027 0.193
## log_ADD_FEE 3.162 2.151 0.097 1.470 0.165 -1.486 7.810
## R_B_COST -0.004 0.002 -0.126 -1.543 0.147 -0.009 0.001
## PH_D -0.573 0.453 -0.147 -1.264 0.228 -1.552 0.406
## SF_RATIO -1.413 0.564 -0.215 -2.507 0.026 -2.631 -0.195
## GRADUAT 0.344 0.238 0.201 1.445 0.172 -0.171 0.860
## OUT_STAT 0.003 0.001 0.286 2.113 0.055 0.000 0.005
## -------------------------------------------------------------------------------------------Здесь можно сравнить по R^2, например.
model_pr_new<-lm(NEW10 ~ AVRCOMB + SF_RATIO, data=pr_data)|>lm.beta()
summary(model_pr_new)##
## Call:
## lm(formula = NEW10 ~ AVRCOMB + SF_RATIO, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.0779 -6.4604 0.1359 4.5401 17.2921
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -133.57926 NA 20.45055 -6.532 7.66e-07 ***
## AVRCOMB 0.17283 0.84064 0.01434 12.055 6.52e-12 ***
## SF_RATIO -1.18802 -0.17520 0.47288 -2.512 0.0188 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.591 on 25 degrees of freedom
## (25 пропущенных наблюдений удалены)
## Multiple R-squared: 0.9309, Adjusted R-squared: 0.9254
## F-statistic: 168.4 on 2 and 25 DF, p-value: 3.104e-15AIC(model_pr_new)## [1] 197.79744.2.8 Построение и анализ графиков
shapiro.test(residuals(model_pr_new))##
## Shapiro-Wilk normality test
##
## data: residuals(model_pr_new)
## W = 0.9525, p-value = 0.2288plot(model_pr_new, 1:5)
У 15-го сильное плечо (ещё видно, как линия идёт к ней)
4.2.9 Поиск outliers по Residuals vs. Deleted Residuals
plot(residuals(model_pr_new), rstudent(model_pr_new))
Все значения не так сильно отличаются, так что посмотрим на outliers по Куку и Махаланобису
4.2.10 Поиск outliers по Cook и Mahalanobis
4.2.10.1 Cook
# Рассчитаем Cook's distance для модели с логарифмами
cooks_vals <- cooks.distance(model_pr_new)
cook_df <- data.frame(
obs = seq_along(cooks_vals),
cook = cooks_vals
) %>%
arrange(desc(cook))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(cook_df))
plot(
x, cook_df$cook,
type = "h",
main = "Cook's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Cook's distance",
xaxt = "n" # убираем автоматическую ось X
)
axis(side = 1, at = x, labels = cook_df$obs, las = 2, cex.axis = 0.8)
abline(h = 1, col = "red", lty = 2)
4.2.10.2 Mahalanobis
# 7.2. Расстояние Махаланобиса (Mahalanobis Distance)
# Расстояние Махаланобиса в пространстве предикторов (регрессоров)
mahalanobis_distance <- mahalanobis(model.matrix(model_pr_new)[,-1],
center = colMeans(model.matrix(model_pr_new)[,-1]),
cov = cov(model.matrix(model_pr_new)[,-1]))
mah_df <- data.frame(
obs = seq_along(mahalanobis_distance),
mah = mahalanobis_distance
) %>%
arrange(desc(mah))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(mah_df))
alpha<-0.05
# Scatter plot для расстояния Махаланобиса
plot(
x, mah_df$mah,
type = "h",
main = "Mahalanobis's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Mahalanobis's distance",
xaxt = "n"
)
axis(side = 1, at = x, labels = mah_df$obs, las = 2, cex.axis = 0.8)
abline(h = qchisq(p = 1 - alpha, df = nrow(mah_df)), col = "red", lty = 2)
Не нашёл outliers (в случае AVRCOMB)
4.2.11 Итог: прогнозирование
set.seed(456) # Другое зерно для генерации новых значений
max_AVRCOMB<-max(pr_data$AVRCOMB, na.rm=TRUE)
max_SF_RATIO<-max(pr_data$SF_RATIO, na.rm=TRUE)
new_data <- data.frame(
AVRCOMB = max_AVRCOMB + 1:8,
SF_RATIO = max_SF_RATIO + 1:8
)
# Прогнозируем значение NEW10
predicted_values <- predict(model_pr_new, newdata = new_data)
cat("Прогнозируемые значения NEW10:\n")## Прогнозируемые значения NEW10:print(predicted_values)## 1 2 3 4 5 6 7 8
## 84.73871 83.72351 82.70832 81.69312 80.67793 79.66273 78.64754 77.63234# 6. Доверительные и предсказательные интервалы
# Доверительный интервал (confidence interval) для среднего значения NEW10
confidence_intervals <- predict(model_pr_new, newdata = new_data, interval = "confidence", level = 0.95)
print("Доверительные интервалы (95%):\n")## [1] "Доверительные интервалы (95%):\n"print(confidence_intervals)## fit lwr upr
## 1 84.73871 68.04653 101.4309
## 2 83.72351 66.08992 101.3571
## 3 82.70832 64.12789 101.2887
## 4 81.69312 62.16121 101.2250
## 5 80.67793 60.19053 101.1653
## 6 79.66273 58.21640 101.1091
## 7 78.64754 56.23924 101.0558
## 8 77.63234 54.25944 101.0052# Предсказательный интервал (prediction interval) для отдельного наблюдения NEW10
prediction_intervals <- predict(model_pr_new, newdata = new_data, interval = "prediction", level = 0.95)
print("Предсказательные интервалы (95%):\n")## [1] "Предсказательные интервалы (95%):\n"print(prediction_intervals)## fit lwr upr
## 1 84.73871 61.86841 107.6090
## 2 83.72351 60.15733 107.2897
## 3 82.70832 58.42552 106.9911
## 4 81.69312 56.67478 106.7115
## 5 80.67793 54.90672 106.4491
## 6 79.66273 53.12282 106.2026
## 7 78.64754 51.32441 105.9707
## 8 77.63234 49.51270 105.7520# Создаем data.frame для прогнозов и интервалов
predictions_df <- data.frame(
AVRCOMB = new_data$AVRCOMB,
SF_RATIO = new_data$SF_RATIO,
Predicted = predicted_values,
Lower_Conf = confidence_intervals[, "lwr"],
Upper_Conf = confidence_intervals[, "upr"],
Lower_Pred = prediction_intervals[, "lwr"],
Upper_Pred = prediction_intervals[, "upr"]
)
predictions_df# Визуализация (только для демонстрации, нужен более сложный график для нескольких предикторов)
ggplot(predictions_df, aes(x = 1:nrow(predictions_df), y = Predicted)) + # Условная ось x, т.к. несколько предикторов
geom_point() +
geom_errorbar(aes(ymin = Lower_Conf, ymax = Upper_Conf), color = "blue", width = 0.2) + # Доверительные интервалы
geom_errorbar(aes(ymin = Lower_Pred, ymax = Upper_Pred), color = "red", width = 0.2) + # Предсказательные интервалы
labs(title = "Прогнозы NEW10 с доверительными и предсказательными интервалами",
x = "Прогноз (номер)",
y = "NEW10") +
scale_x_continuous(breaks = 1:nrow(predictions_df)) + # Подписи для оси x
theme_bw()
# annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.9, label = "Синие: Доверительные интервалы", color = "blue", hjust = 0) +
# annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.8, label = "Красные: Предсказательные интервалы", color = "red", hjust = 0)4.3 Анализ общественных вузов
pub_data<-col_I_regr_split$public
print(paste0("Число наблюдений: ", dim(pub_data)[1]))## [1] "Число наблюдений: 123"colSums(is.na(pub_data))## ...1 log_ADD_FEE BOOK R_B_COST OUT_STAT NEW10
## 0 33 2 0 1 15
## AVRCOMB SF_RATIO PH_D GRADUAT PPIND
## 43 0 5 4 0Много пропусков в AVRCOMB
plot_numeric_pairs(pub_data, "Общественные вузы")
4.3.1 Выбросы “на глаз”
outlier_pub_out_eye<-pub_data[pub_data$OUT_STAT>14001,]
outlier_pub_out_eyeplot_numeric_pairs(pub_data, "Общественные вузы", highlight = list(red=outlier_pub_out_eye))
4.3.2 Построение линейной регрессии
Здесь условия из .tsk файла
Стоит помнить, что у AVRCOMB много NA. В такой модели не будет выбросов “на глаз”
model_pub_new <- lm(NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT,
data = pub_data)|>lm.beta()4.3.3 Интерпретация результатов
summary(model_pub_new)##
## Call:
## lm(formula = NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST +
## PH_D + SF_RATIO + GRADUAT + OUT_STAT, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.462 -7.954 -1.811 5.800 33.454
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -1.957e+02 NA 5.205e+01 -3.760 0.000531 ***
## AVRCOMB 1.124e-01 3.828e-01 4.176e-02 2.691 0.010254 *
## BOOK 3.210e-02 1.521e-01 2.409e-02 1.333 0.190033
## log_ADD_FEE 4.645e+00 2.355e-01 2.416e+00 1.922 0.061536 .
## R_B_COST 4.859e-03 1.906e-01 3.400e-03 1.429 0.160534
## PH_D 4.937e-01 1.328e-01 4.894e-01 1.009 0.318999
## SF_RATIO 7.921e-02 1.507e-02 6.456e-01 0.123 0.902959
## GRADUAT 1.657e-01 1.096e-01 2.419e-01 0.685 0.497348
## OUT_STAT -8.391e-04 -9.592e-02 1.337e-03 -0.628 0.533709
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.48 on 41 degrees of freedom
## (73 пропущенных наблюдений удалены)
## Multiple R-squared: 0.5774, Adjusted R-squared: 0.4949
## F-statistic: 7.002 on 8 and 41 DF, p-value: 8.678e-064.3.4 Построение доверительных эллипсоидов
plot_confidence_ellipse(model_pub_new, c(3, 2), level=0.95)
На примере этой пары: AVRCOMB и SF_RATIO значимы
4.3.5 Множественная корреляция
regressors <- setdiff(names(pub_data), c("...1", "PPIND", "NEW10"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pub_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pub_data[,regressors_full])
)
print(redundancy_table)## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.6083961 0.2875337
## BOOK BOOK 0.4915332 0.2037467
## R_B_COST R_B_COST 0.6691166 0.2178369
## OUT_STAT OUT_STAT 0.7502290 -0.0975559
## AVRCOMB AVRCOMB 0.7530387 0.3874765
## SF_RATIO SF_RATIO 0.5634165 0.0191561
## PH_D PH_D 0.6475179 0.1556248
## GRADUAT GRADUAT 0.7761031 0.1063354Можем убрать OUT_STAT и GRADUAT
AIC(model_pub_new)## [1] 425.9309model_pub1 <- lm(NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO,
data = pub_data)|>lm.beta()
summary(model_pub1)##
## Call:
## lm(formula = NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST +
## PH_D + SF_RATIO, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.863 -7.446 -2.303 7.075 32.472
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -2.035e+02 NA 4.107e+01 -4.955 1.07e-05 ***
## AVRCOMB 1.151e-01 4.081e-01 3.088e-02 3.727 0.000539 ***
## BOOK 3.398e-02 1.619e-01 2.134e-02 1.592 0.118345
## log_ADD_FEE 5.201e+00 2.632e-01 2.298e+00 2.263 0.028499 *
## R_B_COST 3.886e-03 1.546e-01 2.844e-03 1.367 0.178554
## PH_D 5.734e-01 1.550e-01 4.594e-01 1.248 0.218431
## SF_RATIO 1.030e-01 1.985e-02 5.417e-01 0.190 0.849989
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.14 on 45 degrees of freedom
## (71 пропущенное наблюдение удалено)
## Multiple R-squared: 0.5601, Adjusted R-squared: 0.5014
## F-statistic: 9.548 on 6 and 45 DF, p-value: 9.133e-07AIC(model_pub1)## [1] 438.6636AIC увеличился, а значит предпочтительнее выбирать model_pub_new.
regr_high <- setdiff(names(pub_data), c("...1", "PPIND", "NEW10", "OUT_STAT", "GRADUAT"))
regr_fullhigh<-c(regr_high, "NEW10")
redund_high <- data.frame(
Регрессор = regr_high,
Множественная_корреляция = sapply(regr_high, calc_multicollinearity, pub_data[, regr_fullhigh]),
Частная_корреляция = sapply(regr_high, calc_partial_correlation, pub_data[,regr_fullhigh])
)
print(redund_high)## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.5924855 0.31966785
## BOOK BOOK 0.3246989 0.23093274
## R_B_COST R_B_COST 0.5163902 0.19961448
## AVRCOMB AVRCOMB 0.6139954 0.48569728
## SF_RATIO SF_RATIO 0.3206426 0.02834558
## PH_D PH_D 0.6222610 0.182923984.3.6 OLS. Пошаговый AIC.
library(olsrr)
backward_aic_pub<-ols_step_backward_aic(model_pub_new, progress=TRUE, details=TRUE)## Backward Elimination Method
## ---------------------------
##
## Candidate Terms:
##
## 1. AVRCOMB
## 2. BOOK
## 3. log_ADD_FEE
## 4. R_B_COST
## 5. PH_D
## 6. SF_RATIO
## 7. GRADUAT
## 8. OUT_STAT
##
##
## Step => 0
## Model => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 425.9309
##
## Initiating stepwise selection...
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## SF_RATIO 1 423.949 441.157 285.464 0.57725 0.50679
## OUT_STAT 1 424.409 441.617 285.755 0.57334 0.50223
## GRADUAT 1 424.499 441.708 285.813 0.57257 0.50133
## PH_D 1 425.157 442.365 286.230 0.56691 0.49473
## BOOK 1 426.051 443.259 286.802 0.55910 0.48562
## R_B_COST 1 426.362 443.570 287.002 0.55635 0.48241
## log_ADD_FEE 1 428.246 445.454 288.223 0.53932 0.46253
## AVRCOMB 1 432.065 449.273 290.751 0.50275 0.41987
## ------------------------------------------------------------------------
##
## Step => 1
## Removed => SF_RATIO
## Model => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + GRADUAT + OUT_STAT
## AIC => 423.9492
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## GRADUAT 1 422.552 437.849 283.437 0.57212 0.51241
## OUT_STAT 1 422.682 437.978 283.525 0.57101 0.51115
## PH_D 1 423.231 438.527 283.901 0.56627 0.50575
## BOOK 1 424.053 439.349 284.466 0.55908 0.49755
## R_B_COST 1 424.548 439.844 284.807 0.55470 0.49256
## log_ADD_FEE 1 426.379 441.676 286.079 0.53808 0.47363
## AVRCOMB 1 430.065 445.361 288.678 0.50275 0.43336
## ------------------------------------------------------------------------
##
## Step => 2
## Removed => GRADUAT
## Model => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + OUT_STAT
## AIC => 422.5524
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## OUT_STAT 1 420.913 434.297 281.292 0.56902 0.52004
## PH_D 1 421.997 435.381 282.087 0.55957 0.50952
## BOOK 1 422.446 435.830 282.417 0.55560 0.50510
## R_B_COST 1 423.127 436.511 282.920 0.54951 0.49831
## log_ADD_FEE 1 425.924 439.308 284.998 0.52358 0.46945
## AVRCOMB 1 435.531 448.915 292.326 0.42266 0.35706
## ------------------------------------------------------------------------
##
## Step => 3
## Removed => OUT_STAT
## Model => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D
## AIC => 420.9132
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## PH_D 1 420.513 431.985 280.150 0.55501 0.51545
## R_B_COST 1 421.160 432.632 280.657 0.54921 0.50914
## BOOK 1 421.831 433.303 281.185 0.54312 0.50250
## log_ADD_FEE 1 424.596 436.068 283.366 0.51715 0.47423
## AVRCOMB 1 433.721 445.193 290.696 0.42046 0.36895
## ------------------------------------------------------------------------
##
## Step => 4
## Removed => PH_D
## Model => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST
## AIC => 420.5128
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## R_B_COST 1 421.792 431.352 280.561 0.52484 0.49386
## BOOK 1 422.106 431.666 280.822 0.52185 0.49067
## log_ADD_FEE 1 425.733 435.293 283.858 0.48588 0.45235
## AVRCOMB 1 436.611 446.171 293.093 0.36093 0.31925
## ------------------------------------------------------------------------
##
##
## No more variables to be removed.
##
## Variables Removed:
##
## => SF_RATIO
## => GRADUAT
## => OUT_STAT
## => PH_Dforward_aic_pub<-ols_step_forward_aic(model_pub_new, progress=TRUE, details=TRUE)## Forward Selection Method
## ------------------------
##
## Candidate Terms:
##
## 1. AVRCOMB
## 2. BOOK
## 3. log_ADD_FEE
## 4. R_B_COST
## 5. PH_D
## 6. SF_RATIO
## 7. GRADUAT
## 8. OUT_STAT
##
##
## Step => 0
## Model => NEW10 ~ 1
## AIC => 452.9976
##
## Initiating stepwise selection...
##
## Table: Adding New Variables
## -------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## AVRCOMB 1 432.770 438.506 290.015 0.35889 0.34553
## PH_D 1 440.325 446.061 297.005 0.25432 0.23878
## log_ADD_FEE 1 440.669 446.405 297.324 0.24917 0.23353
## GRADUAT 1 440.775 446.511 297.423 0.24757 0.23189
## R_B_COST 1 445.964 451.700 302.247 0.16529 0.14790
## BOOK 1 450.887 456.623 306.843 0.07893 0.05974
## OUT_STAT 1 454.461 460.197 310.192 0.01067 -0.00994
## SF_RATIO 1 454.634 460.370 310.354 0.00724 -0.01344
## -------------------------------------------------------------------------
##
## Step => 1
## Added => AVRCOMB
## Model => NEW10 ~ AVRCOMB
## AIC => 432.7703
##
## Table: Adding New Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## log_ADD_FEE 1 423.754 431.402 281.872 0.48566 0.46377
## R_B_COST 1 426.389 434.038 284.194 0.45782 0.43475
## PH_D 1 427.487 435.135 285.163 0.44579 0.42221
## BOOK 1 431.895 439.543 289.066 0.39471 0.36895
## GRADUAT 1 433.164 440.812 290.193 0.37916 0.35274
## OUT_STAT 1 434.743 442.391 291.599 0.35923 0.33197
## SF_RATIO 1 434.747 442.395 291.602 0.35919 0.33192
## ------------------------------------------------------------------------
##
## Step => 2
## Added => log_ADD_FEE
## Model => NEW10 ~ AVRCOMB + log_ADD_FEE
## AIC => 423.7541
##
## Table: Adding New Variables
## ----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ----------------------------------------------------------------------
## BOOK 1 421.792 431.352 280.561 0.52484 0.49386
## R_B_COST 1 422.106 431.666 280.822 0.52185 0.49067
## PH_D 1 422.181 431.741 280.885 0.52113 0.48990
## OUT_STAT 1 425.569 435.129 283.721 0.48756 0.45414
## GRADUAT 1 425.579 435.140 283.729 0.48745 0.45403
## SF_RATIO 1 425.746 435.306 283.869 0.48574 0.45221
## ----------------------------------------------------------------------
##
## Step => 3
## Added => BOOK
## Model => NEW10 ~ AVRCOMB + log_ADD_FEE + BOOK
## AIC => 421.7919
##
## Table: Adding New Variables
## ----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ----------------------------------------------------------------------
## R_B_COST 1 420.513 431.985 280.150 0.55501 0.51545
## PH_D 1 421.160 432.632 280.657 0.54921 0.50914
## GRADUAT 1 422.984 434.456 282.092 0.53246 0.49091
## OUT_STAT 1 423.748 435.220 282.695 0.52526 0.48306
## SF_RATIO 1 423.783 435.255 282.723 0.52493 0.48270
## ----------------------------------------------------------------------
##
## Step => 4
## Added => R_B_COST
## Model => NEW10 ~ AVRCOMB + log_ADD_FEE + BOOK + R_B_COST
## AIC => 420.5128
##
## Table: Adding New Variables
## ----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ----------------------------------------------------------------------
## PH_D 1 420.913 434.297 281.292 0.56902 0.52004
## OUT_STAT 1 421.997 435.381 282.087 0.55957 0.50952
## GRADUAT 1 422.228 435.613 282.257 0.55753 0.50725
## SF_RATIO 1 422.498 435.882 282.455 0.55514 0.50459
## ----------------------------------------------------------------------
##
##
## No more variables to be added.
##
## Variables Selected:
##
## => AVRCOMB
## => log_ADD_FEE
## => BOOK
## => R_B_COSTПолучили совпадение результатов пошагового forward и пошагового backward.
4.3.6.1 Forward
forward_aic_pub##
##
## Stepwise Summary
## --------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## --------------------------------------------------------------------------
## 0 Base Model 452.998 456.822 309.665 0.00000 0.00000
## 1 AVRCOMB 432.770 438.506 290.015 0.35889 0.34553
## 2 log_ADD_FEE 423.754 431.402 281.872 0.48566 0.46377
## 3 BOOK 421.792 431.352 280.561 0.52484 0.49386
## 4 R_B_COST 420.513 431.985 280.150 0.55501 0.51545
## --------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ----------------------------------------------------------------
## R 0.745 RMSE 14.385
## R-Squared 0.555 MSE 206.930
## Adj. R-Squared 0.515 Coef. Var 42.498
## Pred R-Squared 0.424 AIC 420.513
## MAE 11.265 SBC 431.985
## ----------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## ---------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ---------------------------------------------------------------------
## Regression 12904.403 4 3226.101 14.031 0.0000
## Residual 10346.477 45 229.922
## Total 23250.880 49
## ---------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------------
## (Intercept) -184.284 31.816 -5.792 0.000 -248.364 -120.203
## AVRCOMB 0.136 0.031 0.462 4.430 0.000 0.074 0.197
## log_ADD_FEE 5.845 2.211 0.296 2.644 0.011 1.392 10.297
## BOOK 0.039 0.021 0.185 1.831 0.074 -0.004 0.082
## R_B_COST 0.005 0.003 0.192 1.746 0.088 -0.001 0.011
## ----------------------------------------------------------------------------------------------4.3.6.2 Backward
backward_aic_pub##
##
## Stepwise Summary
## -------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## 0 Full Model 425.931 445.051 287.892 0.57740 0.49494
## 1 SF_RATIO 423.949 441.157 285.030 0.57725 0.50679
## 2 GRADUAT 422.552 437.849 282.885 0.57212 0.51241
## 3 OUT_STAT 420.913 434.297 280.618 0.56902 0.52004
## 4 PH_D 420.513 431.985 279.705 0.55501 0.51545
## -------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ----------------------------------------------------------------
## R 0.745 RMSE 14.385
## R-Squared 0.555 MSE 206.930
## Adj. R-Squared 0.515 Coef. Var 42.498
## Pred R-Squared 0.424 AIC 420.513
## MAE 11.265 SBC 431.985
## ----------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## ---------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ---------------------------------------------------------------------
## Regression 12904.403 4 3226.101 14.031 0.0000
## Residual 10346.477 45 229.922
## Total 23250.880 49
## ---------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------------
## (Intercept) -184.284 31.816 -5.792 0.000 -248.364 -120.203
## AVRCOMB 0.136 0.031 0.462 4.430 0.000 0.074 0.197
## BOOK 0.039 0.021 0.185 1.831 0.074 -0.004 0.082
## log_ADD_FEE 5.845 2.211 0.296 2.644 0.011 1.392 10.297
## R_B_COST 0.005 0.003 0.192 1.746 0.088 -0.001 0.011
## ----------------------------------------------------------------------------------------------model_pub_new<-lm(NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST, data=pub_data)|>lm.beta()
summary(model_pub_new)##
## Call:
## lm(formula = NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST,
## data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.246 -9.548 -1.513 9.069 29.958
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -1.895e+02 NA 2.957e+01 -6.410 5.02e-08 ***
## AVRCOMB 1.261e-01 4.023e-01 2.927e-02 4.308 7.70e-05 ***
## BOOK 4.740e-02 2.089e-01 2.025e-02 2.341 0.0233 *
## log_ADD_FEE 6.998e+00 3.365e-01 2.135e+00 3.278 0.0019 **
## R_B_COST 5.708e-03 2.171e-01 2.666e-03 2.141 0.0371 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.23 on 50 degrees of freedom
## (68 пропущенных наблюдений удалены)
## Multiple R-squared: 0.6167, Adjusted R-squared: 0.5861
## F-statistic: 20.11 on 4 and 50 DF, p-value: 6.357e-10AIC(model_pub_new)## [1] 462.40164.3.7 Построение и анализ графиков
shapiro.test(residuals(model_pub_new))##
## Shapiro-Wilk normality test
##
## data: residuals(model_pub_new)
## W = 0.97633, p-value = 0.3475plot(model_pub_new)
На первом графике Residuals vs. Fitted заметен сильный изгиб (матожидание ошибки ненулевое)
4.3.8 Поиск outliers по Residuals vs. Deleted Residuals
plot(residuals(model_pub_new), rstudent(model_pub_new))
Все значения не так сильно отличаются, так что посмотрим на outliers по Куку и Махаланобису
4.3.9 Поиск outliers по Cook и Mahalanobis
4.3.9.1 Cook
# Рассчитаем Cook's distance для модели с логарифмами
cooks_vals <- cooks.distance(model_pub_new)
cook_df <- data.frame(
obs = seq_along(cooks_vals),
cook = cooks_vals
) %>%
arrange(desc(cook))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(cook_df))
plot(
x, cook_df$cook,
type = "h",
main = "Cook's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Cook's distance",
xaxt = "n" # убираем автоматическую ось X
)
axis(side = 1, at = x, labels = cook_df$obs, las = 2, cex.axis = 0.8)
abline(h = 1, col = "red", lty = 2)
4.3.9.2 Mahalanobis
# 7.2. Расстояние Махаланобиса (Mahalanobis Distance)
# Расстояние Махаланобиса в пространстве предикторов (регрессоров)
mahalanobis_distance <- mahalanobis(model.matrix(model_pub_new)[,-1],
center = colMeans(model.matrix(model_pub_new)[,-1]),
cov = cov(model.matrix(model_pub_new)[,-1]))
mah_df <- data.frame(
obs = seq_along(mahalanobis_distance),
mah = mahalanobis_distance
) %>%
arrange(desc(mah))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(mah_df))
alpha<-0.05
# Scatter plot для расстояния Махаланобиса
plot(
x, mah_df$mah,
type = "h",
main = "Mahalanobis's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Mahalanobis's distance",
xaxt = "n"
)
axis(side = 1, at = x, labels = mah_df$obs, las = 2, cex.axis = 0.8)
abline(h = qchisq(p = 1 - alpha, df = nrow(mah_df)), col = "red", lty = 2)
Не нашёл outliers (в случае AVRCOMB)
4.3.10 Итог: прогнозирование
set.seed(456) # Другое зерно для генерации новых значений
max_AVRCOMB<-max(pub_data$AVRCOMB, na.rm=TRUE)
max_BOOK<-max(pub_data$BOOK, na.rm=TRUE)
max_log_ADD_FEE<-max(pub_data$log_ADD_FEE, na.rm=TRUE)
max_R_B_COST<-max(pub_data$R_B_COST, na.rm=TRUE)
new_data <- data.frame(
AVRCOMB = max_AVRCOMB - seq(50, 120, by = 10),
BOOK = max_BOOK - seq(10, 80, by = 10),
log_ADD_FEE = max_log_ADD_FEE - seq(0.01, 0.08, by=0.01),
R_B_COST = max_R_B_COST - seq(10, 80, by=10)
)
# Прогнозируем значение NEW10
predicted_values <- predict(model_pub_new, newdata = new_data)
cat("Прогнозируемые значения NEW10:\n")## Прогнозируемые значения NEW10:print(predicted_values)## 1 2 3 4 5 6 7 8
## 96.97503 95.11299 93.25095 91.38891 89.52687 87.66483 85.80280 83.94076# 6. Доверительные и предсказательные интервалы
# Доверительный интервал (confidence interval) для среднего значения NEW10
confidence_intervals <- predict(model_pub_new, newdata = new_data, interval = "confidence", level = 0.95)
print("Доверительные интервалы (95%):\n")## [1] "Доверительные интервалы (95%):\n"print(confidence_intervals)## fit lwr upr
## 1 96.97503 81.95568 111.99437
## 2 95.11299 80.48533 109.74065
## 3 93.25095 78.99581 107.50609
## 4 91.38891 77.48558 105.29224
## 5 89.52687 75.95304 103.10071
## 6 87.66483 74.39651 100.93315
## 7 85.80280 72.81432 98.79127
## 8 83.94076 71.20476 96.67675# Предсказательный интервал (prediction interval) для отдельного наблюдения NEW10
prediction_intervals <- predict(model_pub_new, newdata = new_data, interval = "prediction", level = 0.95)
print("Предсказательные интервалы (95%):\n")## [1] "Предсказательные интервалы (95%):\n"print(prediction_intervals)## fit lwr upr
## 1 96.97503 62.89617 131.0539
## 2 95.11299 61.20493 129.0210
## 3 93.25095 59.50192 127.0000
## 4 91.38891 57.78697 124.9909
## 5 89.52687 56.05992 122.9938
## 6 87.66483 54.32063 121.0090
## 7 85.80280 52.56895 119.0366
## 8 83.94076 50.80477 117.0767predictions_df <- data.frame(
AVRCOMB = new_data$AVRCOMB,
BOOK = new_data$BOOK,
log_ADD_FEE = new_data$log_ADD_FEE,
R_B_COST = new_data$R_B_COST,
Predicted = predicted_values,
Lower_Conf = confidence_intervals[, "lwr"],
Upper_Conf = confidence_intervals[, "upr"],
Lower_Pred = prediction_intervals[, "lwr"],
Upper_Pred = prediction_intervals[, "upr"]
)
predictions_dfggplot(predictions_df, aes(x = 1:nrow(predictions_df), y = Predicted)) +
geom_point() +
geom_errorbar(aes(ymin = Lower_Conf, ymax = Upper_Conf), color = "blue", width = 0.2) + # Доверительные интервалы
geom_errorbar(aes(ymin = Lower_Pred, ymax = Upper_Pred), color = "red", width = 0.2) + # Предсказательные интервалы
labs(title = "Прогнозы NEW10 с доверительными и предсказательными интервалами",
x = "Прогноз (номер)",
y = "NEW10") +
scale_x_continuous(breaks = 1:nrow(predictions_df)) + # Подписи для оси x
theme_bw()
4.4 Исправления.
Важные пункты для исправления:
Использование AVRCOMB (особенно в случае частных вузов) не интересно. Поэтому уберём её из рассматриваемых признаков.
Доверительный эллипсоид некорректно строится, поэтому нужно использовать другую функцию построения доверительного эллипсоида.
Рассмотреть прологарифмированный NEW10 для общественных вузов
4.5 Подготовка данных
4.5.1 Выбор признаков
Рассматриваемые наблюдения: 1. PPIND (для группировки) 2. ADD_FEE - дополнительные выплаты 3. BOOK - плата за бронирование 4. NEW10 (зависимая переменная, также NEW25) - студентов из топ 10% своей старшей школы 5. AVRCOMB - средний комбинированный балл (лишний признак) 6. SF_RATIO - Student/Faculty ratio 7. PH_D - Pct. of faculty with Ph.D.’s 8. GRADUAT - Graduation rate 9. R_B_COST - Room and board costs 10. OUT_STAT - Out-of-state tuition
colNames<-c("...1", "ADD_FEE", "BOOK", "R_B_COST","OUT_STAT","NEW10", "SF_RATIO","PH_D","GRADUAT","PPIND")
col_I_regr<-col_I_sn[c(colNames)]
col_I_regr<-col_I_regr|>mutate(ADD_FEE = log(ADD_FEE))|>rename(log_ADD_FEE=ADD_FEE)
col_I_regr_split<-split(col_I_regr, col_I_regr$PPIND)# Автор: -Николай-
library(GGally)
library(ggplot2)
library(dplyr)
library(rlang)
plot_numeric_pairs <- function(df, title = NULL, highlight = NULL) {
# Выбираем имена всех numeric-столбцов
numeric_cols <- df %>% dplyr::select(dplyr::where(is.numeric)) %>% names()
# Функция для нижних панелей с возможностью выделения точек разными цветами
custom_points <- function(data, mapping, ...) {
# Отрисовываем базовые точки синим цветом
p <- ggplot(data, mapping) +
geom_point(alpha = 0.6, size = 2, color = "blue", ...)
if (!is.null(highlight)) {
# Если highlight не является списком, превращаем его в список с одним элементом и цветом "red"
if (!is.list(highlight) || inherits(highlight, "data.frame")) {
highlight <- list(red = highlight)
}
# Получаем имена переменных из mapping с помощью as_label()
xvar <- as_label(mapping$x)
yvar <- as_label(mapping$y)
# Для каждого набора точек в списке highlight
for (col in names(highlight)) {
hl_df <- highlight[[col]]
# Проверяем, что hl_df содержит нужные переменные
if(all(c(xvar, yvar) %in% names(hl_df))) {
# Находим совпадающие строки между данными панели и hl_df
data_highlight <- semi_join(data, hl_df, by = c(xvar, yvar))
# Накладываем выделенные точки нужного цвета
p <- p + geom_point(data = data_highlight, mapping = mapping,
alpha = 0.6, size = 2, color = col, ...)
}
}
}
p
}
# Строим матрицу парных графиков
ggpairs(df,
columns = numeric_cols,
columnLabels = numeric_cols,
title = title,
upper = list(continuous = wrap("cor", size = 3)),
lower = list(continuous = custom_points),
diag = list(continuous = wrap("barDiag", bins = 10, fill = "lightblue", color = "black"))
) +
theme_bw() +
theme(
strip.text.x = element_text(size = 7, angle = 0, hjust = 0),
strip.text.y = element_text(size = 8, angle = 0, hjust = 0),
axis.text = element_text(size = 6),
plot.margin = unit(c(1, 2, 1, 1), "cm")
)
}4.6 Общий парный график
ggpairs(col_I_regr[,-which(names(col_I_regr)=="...1")], aes(alpha = 0.5, color=PPIND),
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("cor", size = 2)))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Комментарий: замечаем, что для публичных вузов NEW10 имеет выраженный правый хвост. Поэтому стоит использовать log_NEW10 для публичных вузов.
При постановке вопроса о построении линейной модели с категориальным признаком стоит обратить внимание на то, близки ли корреляции разных групп. Если они отличаются, то лучше рассматривать каждую категорию по отдельности. Если близки, то вместе. Это связано с тем, что разные корреляции по разному влияют на общую модель, будут возникать неоднородности.
4.7 Анализ частных вузов
pr_data<-col_I_regr_split$private
print(paste0("Число наблюдений: ", dim(pr_data)[1]))## [1] "Число наблюдений: 53"colSums(is.na(pr_data))## ...1 log_ADD_FEE BOOK R_B_COST OUT_STAT NEW10
## 0 7 0 0 0 1
## SF_RATIO PH_D GRADUAT PPIND
## 0 4 1 04.7.1 График
plot_numeric_pairs(pr_data, "Частные вузы")
4.7.2 Выбросы “на глаз”
outlier_pr_book_eye<-pr_data[pr_data$BOOK>875,]
outlier_pr_out_eye<-pr_data[pr_data$OUT_STAT<5000,]
outlier_pr_book_eyeoutlier_pr_out_eyeplot_numeric_pairs(pr_data, "Частные вузы", list(red = outlier_pr_book_eye, green = outlier_pr_out_eye))
4.7.3 Линейная модель
library(lm.beta)
model_pr_full <- lm(NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT,
data = pr_data) |> lm.beta()
summary(model_pr_full)##
## Call:
## lm(formula = NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO +
## GRADUAT + OUT_STAT, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.1636 -6.8538 0.0723 7.2414 25.5012
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -9.289705 NA 38.856630 -0.239 0.81252
## BOOK 0.002300 0.015211 0.013944 0.165 0.86999
## log_ADD_FEE 1.877765 0.063416 2.978846 0.630 0.53280
## R_B_COST -0.002897 -0.120242 0.002971 -0.975 0.33658
## PH_D -0.219550 -0.069638 0.508017 -0.432 0.66843
## SF_RATIO -1.772400 -0.324548 0.604579 -2.932 0.00608 **
## GRADUAT 1.069124 0.642349 0.215538 4.960 2.07e-05 ***
## OUT_STAT 0.001308 0.192473 0.001098 1.192 0.24180
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.4 on 33 degrees of freedom
## (12 пропущенных наблюдений удалены)
## Multiple R-squared: 0.73, Adjusted R-squared: 0.6727
## F-statistic: 12.75 on 7 and 33 DF, p-value: 8.508e-084.7.4 Доверительные эллипсы
plot_confidence_ellipse(model_pr_full, c(2, 3))
plot_confidence_ellipse(model_pr_full, c(4, 5))
plot_confidence_ellipse(model_pr_full, c(6, 7))
plot_confidence_ellipse(model_pr_full, c(8, 2))
# Автор: -Евгений-
# Вычисление множественной корреляции
calc_multicollinearity <- function(var, data) {
X <- as.matrix(data[, setdiff(names(data), var)]) # Матрица остальных регрессоров
y <- as.matrix(data[[var]]) # Текущий регрессор как зависимая переменная
model <- lm(y ~ X)
sqrt(summary(model)$r.squared)
}
calc_partial_correlation <- function(var, data) {
# Остатки регрессии текущего регрессора на остальные регрессоры
data<-data|>na.omit()
model_X <- lm(data[[var]] ~ ., data = data[, !names(data) %in% c(var, "NEW10")])
res_X <- residuals(model_X)
# Остатки регрессии Y на остальные регрессоры
model_Y <- lm(NEW10 ~ ., data = data[, !names(data) %in% var])
res_Y <- residuals(model_Y)
cor(res_X, res_Y, use = "complete.obs") # Корреляция остатков
}
regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)
print(redundancy_table)## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.4485260 0.10907796
## BOOK BOOK 0.1965671 0.02870252
## R_B_COST R_B_COST 0.6906122 -0.16735732
## OUT_STAT OUT_STAT 0.8361740 0.20315778
## SF_RATIO SF_RATIO 0.6858132 -0.45456011
## PH_D PH_D 0.8286465 -0.07501937
## GRADUAT GRADUAT 0.8488216 0.65354784Комментарий: Стоит исключить PH_D из-за высокой множественной и низкой частной корреляций (также можно исключить BOOK из-за низкой частной корреляции).
AIC(model_pr_full)## [1] 338.2841summary(model_pr_full)##
## Call:
## lm(formula = NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO +
## GRADUAT + OUT_STAT, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.1636 -6.8538 0.0723 7.2414 25.5012
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -9.289705 NA 38.856630 -0.239 0.81252
## BOOK 0.002300 0.015211 0.013944 0.165 0.86999
## log_ADD_FEE 1.877765 0.063416 2.978846 0.630 0.53280
## R_B_COST -0.002897 -0.120242 0.002971 -0.975 0.33658
## PH_D -0.219550 -0.069638 0.508017 -0.432 0.66843
## SF_RATIO -1.772400 -0.324548 0.604579 -2.932 0.00608 **
## GRADUAT 1.069124 0.642349 0.215538 4.960 2.07e-05 ***
## OUT_STAT 0.001308 0.192473 0.001098 1.192 0.24180
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.4 on 33 degrees of freedom
## (12 пропущенных наблюдений удалены)
## Multiple R-squared: 0.73, Adjusted R-squared: 0.6727
## F-statistic: 12.75 on 7 and 33 DF, p-value: 8.508e-08model_pr_nono<-lm(NEW10 ~ log_ADD_FEE + R_B_COST + SF_RATIO + GRADUAT + OUT_STAT, data=pr_data)|>lm.beta()
summary(model_pr_nono)##
## Call:
## lm(formula = NEW10 ~ log_ADD_FEE + R_B_COST + SF_RATIO + GRADUAT +
## OUT_STAT, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.1896 -7.8512 0.4492 6.5775 29.4475
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -1.254e+01 NA 2.304e+01 -0.544 0.5895
## log_ADD_FEE 1.313e+00 4.221e-02 2.810e+00 0.467 0.6430
## R_B_COST -3.726e-03 -1.462e-01 2.743e-03 -1.358 0.1823
## SF_RATIO -1.872e+00 -3.275e-01 5.701e-01 -3.284 0.0022 **
## GRADUAT 1.019e+00 5.877e-01 1.880e-01 5.420 3.56e-06 ***
## OUT_STAT 1.228e-03 1.710e-01 9.368e-04 1.311 0.1976
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.49 on 38 degrees of freedom
## (9 пропущенных наблюдений удалены)
## Multiple R-squared: 0.726, Adjusted R-squared: 0.6899
## F-statistic: 20.13 on 5 and 38 DF, p-value: 9.083e-10AIC(model_pr_nono)## [1] 361.386Комментарий: замечаем, что после того, как убрали BOOK и PH_D увеличился Adj.R^2, немного уменьшился R^2, уменьшился p-value, но увеличился AIC.
regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10", "PH_D", "BOOK"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)
print(redundancy_table)## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.3480743 0.07558012
## R_B_COST R_B_COST 0.6376041 -0.21521296
## OUT_STAT OUT_STAT 0.7708944 0.20809285
## SF_RATIO SF_RATIO 0.6597111 -0.47020518
## GRADUAT GRADUAT 0.8087224 0.66027976Попробуем убрать log_ADD_FEE
model_pr1<-lm(NEW10 ~ GRADUAT + R_B_COST + SF_RATIO + OUT_STAT, data=pr_data)|>lm.beta()
summary(model_pr1)##
## Call:
## lm(formula = NEW10 ~ GRADUAT + R_B_COST + SF_RATIO + OUT_STAT,
## data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.041 -9.947 -1.776 8.355 49.615
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 14.1924930 NA 21.7789414 0.652 0.51786
## GRADUAT 0.9698293 0.6097925 0.1999565 4.850 1.45e-05 ***
## R_B_COST -0.0038847 -0.1556376 0.0029142 -1.333 0.18908
## SF_RATIO -1.6934700 -0.3080606 0.6201252 -2.731 0.00893 **
## OUT_STAT 0.0003467 0.0543364 0.0010089 0.344 0.73268
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.36 on 46 degrees of freedom
## (2 пропущенных наблюдений удалены)
## Multiple R-squared: 0.6291, Adjusted R-squared: 0.5968
## F-statistic: 19.5 on 4 and 46 DF, p-value: 1.918e-09AIC(model_pr1)## [1] 430.1013Комментарий: AIC увеличилась, а R^2 сильно уменьшились, p-value увеличился.
regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10", "PH_D", "BOOK", "log_ADD_FEE"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)
print(redundancy_table)## Регрессор Множественная_корреляция Частная_корреляция
## R_B_COST R_B_COST 0.6561089 -0.19285565
## OUT_STAT OUT_STAT 0.8235937 0.05060282
## SF_RATIO SF_RATIO 0.6743579 -0.37350259
## GRADUAT GRADUAT 0.8139325 0.58168914AIC(model_pr_full)## [1] 338.2841AIC(model_pr_nono)## [1] 361.386AIC(model_pr1)## [1] 430.10134.7.5 OLS. Пошаговый AIC.
4.7.5.1 Backward
library(olsrr)
backward_aic<-ols_step_backward_aic(model_pr_full, progress=TRUE, details=TRUE)## Backward Elimination Method
## ---------------------------
##
## Candidate Terms:
##
## 1. BOOK
## 2. log_ADD_FEE
## 3. R_B_COST
## 4. PH_D
## 5. SF_RATIO
## 6. GRADUAT
## 7. OUT_STAT
##
##
## Step => 0
## Model => NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 338.2841
##
## Initiating stepwise selection...
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## BOOK 1 336.318 350.026 223.228 0.72977 0.68209
## PH_D 1 336.516 350.224 223.348 0.72847 0.68055
## log_ADD_FEE 1 336.775 350.483 223.505 0.72674 0.67852
## R_B_COST 1 337.449 351.157 223.916 0.72221 0.67319
## OUT_STAT 1 338.012 351.721 224.262 0.71837 0.66867
## SF_RATIO 1 345.774 359.483 229.220 0.65968 0.59962
## GRADUAT 1 359.125 372.833 238.570 0.52868 0.44551
## ------------------------------------------------------------------------
##
## Step => 1
## Removed => BOOK
## Model => NEW10 ~ log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 336.3179
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## PH_D 1 334.535 346.530 220.890 0.72834 0.68953
## log_ADD_FEE 1 334.821 346.816 221.080 0.72644 0.68736
## R_B_COST 1 335.500 347.495 221.536 0.72187 0.68213
## OUT_STAT 1 336.045 348.040 221.902 0.71815 0.67789
## SF_RATIO 1 343.786 355.781 227.258 0.65958 0.61094
## GRADUAT 1 357.284 369.279 237.270 0.52684 0.45925
## ------------------------------------------------------------------------
##
## Step => 2
## Removed => PH_D
## Model => NEW10 ~ log_ADD_FEE + R_B_COST + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 334.5351
##
## Table: Removing Existing Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## log_ADD_FEE 1 332.891 343.172 218.679 0.72597 0.69552
## R_B_COST 1 333.513 343.794 219.132 0.72178 0.69087
## OUT_STAT 1 334.111 344.392 219.570 0.71770 0.68633
## SF_RATIO 1 341.927 352.208 225.394 0.65840 0.62045
## GRADUAT 1 358.480 368.761 238.448 0.48849 0.43166
## ------------------------------------------------------------------------
##
## Step => 3
## Removed => log_ADD_FEE
## Model => NEW10 ~ R_B_COST + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 332.891
##
## Table: Removing Existing Variables
## ----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ----------------------------------------------------------------------
## R_B_COST 1 331.621 340.189 216.806 0.72105 0.69843
## OUT_STAT 1 332.295 340.863 217.339 0.71642 0.69343
## SF_RATIO 1 339.927 348.495 223.430 0.65840 0.63071
## GRADUAT 1 358.030 366.598 238.458 0.46878 0.42571
## ----------------------------------------------------------------------
##
## Step => 4
## Removed => R_B_COST
## Model => NEW10 ~ SF_RATIO + GRADUAT + OUT_STAT
## AIC => 331.6207
##
## Table: Removing Existing Variables
## ----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ----------------------------------------------------------------------
## OUT_STAT 1 330.381 337.236 215.046 0.71582 0.70087
## SF_RATIO 1 338.121 344.975 221.632 0.65678 0.63872
## GRADUAT 1 358.218 365.073 239.135 0.43966 0.41017
## ----------------------------------------------------------------------
##
## Step => 5
## Removed => OUT_STAT
## Model => NEW10 ~ SF_RATIO + GRADUAT
## AIC => 330.3813
##
## Table: Removing Existing Variables
## ----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ----------------------------------------------------------------------
## SF_RATIO 1 339.980 345.121 223.161 0.62291 0.61324
## GRADUAT 1 365.570 370.711 246.575 0.29611 0.27806
## ----------------------------------------------------------------------
##
##
## No more variables to be removed.
##
## Variables Removed:
##
## => BOOK
## => PH_D
## => log_ADD_FEE
## => R_B_COST
## => OUT_STATbackward_aic##
##
## Stepwise Summary
## --------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## --------------------------------------------------------------------------
## 0 Full Model 338.284 353.706 225.692 0.73000 0.67272
## 1 BOOK 336.318 350.026 222.763 0.72977 0.68209
## 2 PH_D 334.535 346.530 220.181 0.72834 0.68953
## 3 log_ADD_FEE 332.891 343.172 217.888 0.72597 0.69552
## 4 R_B_COST 331.621 340.189 216.109 0.72105 0.69843
## 5 OUT_STAT 330.381 337.236 214.490 0.71582 0.70087
## --------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ----------------------------------------------------------------
## R 0.846 RMSE 12.336
## R-Squared 0.716 MSE 152.183
## Adj. R-Squared 0.701 Coef. Var 23.175
## Pred R-Squared 0.672 AIC 330.381
## MAE 10.318 SBC 337.236
## ----------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## ---------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ---------------------------------------------------------------------
## Regression 15717.001 2 7858.501 47.86 0.0000
## Residual 6239.487 38 164.197
## Total 21956.488 40
## ---------------------------------------------------------------------
##
## Parameter Estimates
## ------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ------------------------------------------------------------------------------------------
## (Intercept) -18.639 14.637 -1.273 0.211 -48.270 10.992
## SF_RATIO -1.760 0.499 -0.322 -3.525 0.001 -2.770 -0.749
## GRADUAT 1.140 0.152 0.685 7.492 0.000 0.832 1.448
## ------------------------------------------------------------------------------------------4.7.5.2 Forward
forward_aic<-ols_step_forward_aic(model_pr_full, progress=TRUE, details=TRUE)## Forward Selection Method
## ------------------------
##
## Candidate Terms:
##
## 1. BOOK
## 2. log_ADD_FEE
## 3. R_B_COST
## 4. PH_D
## 5. SF_RATIO
## 6. GRADUAT
## 7. OUT_STAT
##
##
## Step => 0
## Model => NEW10 ~ 1
## AIC => 377.966
##
## Initiating stepwise selection...
##
## Table: Adding New Variables
## -------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## GRADUAT 1 339.980 345.121 223.161 0.62291 0.61324
## PH_D 1 358.045 363.186 239.617 0.41413 0.39910
## OUT_STAT 1 362.074 367.215 243.334 0.35363 0.33706
## SF_RATIO 1 365.570 370.711 246.575 0.29611 0.27806
## R_B_COST 1 378.009 383.150 258.223 0.04661 0.02216
## BOOK 1 379.594 384.735 259.720 0.00903 -0.01638
## log_ADD_FEE 1 379.787 384.928 259.903 0.00435 -0.02118
## -------------------------------------------------------------------------
##
## Step => 1
## Added => GRADUAT
## Model => NEW10 ~ GRADUAT
## AIC => 339.9804
##
## Table: Adding New Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## SF_RATIO 1 330.381 337.236 215.046 0.71582 0.70087
## OUT_STAT 1 338.121 344.975 221.632 0.65678 0.63872
## PH_D 1 340.187 347.041 223.400 0.63905 0.62005
## R_B_COST 1 341.456 348.310 224.490 0.62770 0.60810
## log_ADD_FEE 1 341.939 348.793 224.905 0.62329 0.60346
## BOOK 1 341.973 348.827 224.934 0.62297 0.60313
## ------------------------------------------------------------------------
##
## Step => 2
## Added => SF_RATIO
## Model => NEW10 ~ GRADUAT + SF_RATIO
## AIC => 330.3813
##
## Table: Adding New Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## OUT_STAT 1 331.621 340.189 216.806 0.72105 0.69843
## log_ADD_FEE 1 332.265 340.833 217.315 0.71663 0.69366
## PH_D 1 332.272 340.840 217.320 0.71658 0.69360
## R_B_COST 1 332.295 340.863 217.339 0.71642 0.69343
## BOOK 1 332.321 340.889 217.359 0.71624 0.69323
## ------------------------------------------------------------------------
##
##
## No more variables to be added.
##
## Variables Selected:
##
## => GRADUAT
## => SF_RATIOforward_aic##
##
## Stepwise Summary
## -------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## 0 Base Model 377.966 381.393 259.401 0.00000 0.00000
## 1 GRADUAT 339.980 345.121 223.161 0.62291 0.61324
## 2 SF_RATIO 330.381 337.236 215.046 0.71582 0.70087
## -------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ----------------------------------------------------------------
## R 0.846 RMSE 12.336
## R-Squared 0.716 MSE 152.183
## Adj. R-Squared 0.701 Coef. Var 23.175
## Pred R-Squared 0.672 AIC 330.381
## MAE 10.318 SBC 337.236
## ----------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## ---------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ---------------------------------------------------------------------
## Regression 15717.001 2 7858.501 47.86 0.0000
## Residual 6239.487 38 164.197
## Total 21956.488 40
## ---------------------------------------------------------------------
##
## Parameter Estimates
## ------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ------------------------------------------------------------------------------------------
## (Intercept) -18.639 14.637 -1.273 0.211 -48.270 10.992
## GRADUAT 1.140 0.152 0.685 7.492 0.000 0.832 1.448
## SF_RATIO -1.760 0.499 -0.322 -3.525 0.001 -2.770 -0.749
## ------------------------------------------------------------------------------------------Комментарий: строим Backward и Forward AIC. Backward AIC исключил всё, кроме SF_RATIO и GRADUAT. Forward AIC включил только GRADUAT и SF_RATIO. Оба пошаговых алгоритма дали один и тот же результат.
Интерпретация параметров: число лучших 10% студентов NEW10 увеличивается с увеличением выпускников GRADUAT и уменьшается с числом студентов на преподавателя SF_RATIO.
Здесь можно сравнить по R^2, например.
4.7.6 Построение и анализ графиков
model_pr_full<-lm(NEW10 ~ GRADUAT + SF_RATIO, data=pr_data)
summary(model_pr_full)##
## Call:
## lm(formula = NEW10 ~ GRADUAT + SF_RATIO, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.763 -8.954 -2.873 9.404 52.149
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.1676 16.4411 -0.314 0.75465
## GRADUAT 0.9740 0.1626 5.990 2.6e-07 ***
## SF_RATIO -1.5182 0.5621 -2.701 0.00952 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.35 on 48 degrees of freedom
## (2 пропущенных наблюдений удалены)
## Multiple R-squared: 0.6133, Adjusted R-squared: 0.5972
## F-statistic: 38.07 on 2 and 48 DF, p-value: 1.247e-10plot(model_pr_full, 1:5)
У 51-го сильное плечо (ещё видно, как линия идёт к ней)
4.7.7 Поиск outliers по Residuals vs. Deleted Residuals
plot(residuals(model_pr_full), rstudent(model_pr_full))
Комментарий: замечаем точку (51), которая сильно влияет на линейную модель.
influence.measures(model_pr_full)## Influence measures of
## lm(formula = NEW10 ~ GRADUAT + SF_RATIO, data = pr_data) :
##
## dfb.1_ dfb.GRAD dfb.SF_R dffit cov.r cook.d hat inf
## 1 0.156252 -0.045588 -2.35e-01 0.4007 0.816 4.95e-02 0.0317
## 2 -0.142405 0.171118 7.76e-02 0.2199 1.057 1.61e-02 0.0498
## 3 0.007595 -0.007708 -9.15e-05 0.0145 1.099 7.16e-05 0.0308
## 4 0.007442 0.011854 -7.56e-02 -0.1148 1.115 4.46e-03 0.0579
## 5 -0.062742 0.118154 -4.21e-02 0.2030 1.079 1.38e-02 0.0554
## 6 0.209642 -0.134914 -3.80e-01 -0.4260 1.077 5.98e-02 0.1036
## 7 -0.010436 0.018846 -3.62e-02 -0.0813 1.093 2.24e-03 0.0360
## 8 -0.156587 0.114409 8.70e-02 -0.3109 0.849 3.03e-02 0.0231
## 9 -0.157795 0.122821 1.25e-01 -0.2056 1.037 1.41e-02 0.0387
## 10 0.025289 -0.043555 6.54e-03 -0.0741 1.102 1.86e-03 0.0411
## 12 -0.031867 0.027892 2.43e-02 -0.0335 1.206 3.81e-04 0.1170 *
## 13 0.007036 -0.029448 3.02e-02 -0.0739 1.117 1.85e-03 0.0523
## 14 0.148220 -0.160722 -2.28e-02 0.2043 1.101 1.40e-02 0.0676
## 15 -0.440485 0.369605 3.53e-01 -0.4841 0.980 7.54e-02 0.0802
## 16 -0.020974 0.058786 -3.49e-02 0.1383 1.069 6.44e-03 0.0357
## 17 -0.013268 -0.014281 4.90e-02 -0.0892 1.099 2.70e-03 0.0421
## 18 -0.165762 0.166156 1.42e-01 0.2044 1.139 1.41e-02 0.0899
## 19 -0.037008 0.028235 2.21e-02 -0.0639 1.082 1.39e-03 0.0251
## 20 -0.014143 0.000378 3.18e-02 -0.0433 1.136 6.37e-04 0.0640
## 21 -0.169351 0.147972 1.93e-01 0.2218 1.183 1.66e-02 0.1200
## 22 0.024521 -0.044146 -1.97e-02 -0.1356 1.034 6.15e-03 0.0219
## 24 -0.233928 0.276537 1.44e-01 0.3683 0.934 4.35e-02 0.0451
## 25 -0.081062 0.088996 1.72e-02 -0.1024 1.203 3.56e-03 0.1191 *
## 26 0.071188 -0.084065 -4.77e-02 -0.1185 1.086 4.74e-03 0.0399
## 27 -0.182509 0.145727 1.67e-01 -0.2044 1.152 1.41e-02 0.0979
## 28 -0.067432 -0.008354 1.65e-01 -0.2386 1.063 1.90e-02 0.0563
## 29 -0.008973 0.028263 -2.49e-02 0.0639 1.127 1.39e-03 0.0585
## 30 -0.110949 0.144068 3.66e-02 0.1853 1.086 1.15e-02 0.0551
## 31 0.084879 -0.059446 -1.54e-01 -0.1872 1.099 1.18e-02 0.0624
## 32 0.022119 -0.040726 1.33e-02 -0.0689 1.121 1.61e-03 0.0550
## 33 -0.016792 0.037689 -1.59e-02 0.0784 1.095 2.09e-03 0.0372
## 34 0.154239 -0.125010 -2.45e-01 -0.3179 0.989 3.30e-02 0.0484
## 35 -0.052304 0.065696 -2.44e-02 -0.1143 1.110 4.43e-03 0.0545
## 36 0.189408 -0.151498 -1.43e-01 0.2417 1.014 1.93e-02 0.0395
## 37 -0.029891 -0.001702 1.42e-01 0.2076 1.068 1.44e-02 0.0514
## 38 -0.105057 0.106146 9.90e-03 -0.1824 1.035 1.11e-02 0.0329
## 39 -0.060391 0.032813 7.43e-02 -0.1045 1.091 3.70e-03 0.0397
## 40 -0.018011 0.057624 -5.00e-02 0.1320 1.105 5.89e-03 0.0543
## 41 0.373278 -0.291574 -3.66e-01 0.4214 1.143 5.90e-02 0.1335
## 42 0.041672 -0.031904 -2.39e-02 0.0735 1.077 1.83e-03 0.0247
## 43 -0.091131 0.106561 -1.98e-02 -0.1722 1.081 9.97e-03 0.0494
## 44 0.288062 -0.292022 -2.68e-01 -0.4058 0.948 5.29e-02 0.0556
## 45 -0.015808 0.050662 -3.76e-02 0.1224 1.083 5.06e-03 0.0394
## 46 -0.062436 0.092305 4.33e-03 0.1352 1.091 6.17e-03 0.0466
## 47 0.024276 -0.009007 -3.60e-02 0.0538 1.102 9.82e-04 0.0377
## 48 -0.073400 0.113328 1.46e-02 0.1985 1.019 1.31e-02 0.0315
## 49 0.000713 0.003235 -2.22e-02 -0.0423 1.099 6.09e-04 0.0335
## 50 -0.004980 0.008760 -1.19e-02 -0.0260 1.121 2.29e-04 0.0506
## 51 1.078464 -1.600184 8.12e-01 2.5174 0.450 1.48e+00 0.2292 *
## 52 0.047770 -0.034990 -1.09e-01 -0.1676 1.051 9.41e-03 0.0348
## 53 0.027573 -0.015210 -5.76e-02 -0.0660 1.184 1.48e-03 0.1022row.names(model_pr_full$model)## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "12" "13" "14" "15" "16"
## [16] "17" "18" "19" "20" "21" "22" "24" "25" "26" "27" "28" "29" "30" "31" "32"
## [31] "33" "34" "35" "36" "37" "38" "39" "40" "41" "42" "43" "44" "45" "46" "47"
## [46] "48" "49" "50" "51" "52" "53"outlier_cook<-pr_data|>dplyr::filter(NEW10==48&SF_RATIO==20.5)outlier_cookКомментарий: выбросом оказался Brigham Young University.
4.7.8 Поиск outliers по Cook и Mahalanobis
4.7.8.1 Cook
# Рассчитаем Cook's distance для модели с логарифмами
cooks_vals <- cooks.distance(model_pr_full)
cook_df <- data.frame(
obs = row.names(model_pr_full$model),
cook = cooks_vals
) %>%
arrange(desc(cook))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(cook_df))
plot(
x, cook_df$cook,
type = "h",
main = "Cook's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Cook's distance",
xaxt = "n" # убираем автоматическую ось X
)
axis(side = 1, at = x, labels = cook_df$obs, las = 2, cex.axis = 0.8)
abline(h = 1, col = "red", lty = 2)
Комментарий: высокие значения у 29 и 38.
4.7.8.2 Mahalanobis
# 7.2. Расстояние Махаланобиса (Mahalanobis Distance)
# Расстояние Махаланобиса в пространстве предикторов (регрессоров)
mahalanobis_distance <- mahalanobis(model.matrix(model_pr_full)[,-1],
center = colMeans(model.matrix(model_pr_full)[,-1]),
cov = cov(model.matrix(model_pr_full)[,-1]))
mah_df <- data.frame(
obs = row.names(model_pr_full$model),
mah = mahalanobis_distance
) %>%
arrange(desc(mah))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(mah_df))
alpha<-0.05
# Scatter plot для расстояния Махаланобиса
plot(
x, mah_df$mah,
type = "h",
main = "Mahalanobis's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Mahalanobis's distance",
xaxt = "n"
)
axis(side = 1, at = x, labels = mah_df$obs, las = 2, cex.axis = 0.8)
abline(h = qchisq(p = 1 - alpha, df = nrow(mah_df)), col = "red", lty = 2)
Комментарий: 29-ое наблюдение имеет большое расстояние по Махаланобису.
4.7.9 Линейная модель без выброса
pr_data_outlier<-pr_data[-c(51),]plot_numeric_pairs(pr_data_outlier)
model_pr_outlier<-lm(NEW10~GRADUAT + SF_RATIO, data=pr_data_outlier)
summary(model_pr_full|>lm.beta())##
## Call:
## lm(formula = NEW10 ~ GRADUAT + SF_RATIO, data = pr_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.763 -8.954 -2.873 9.404 52.149
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -5.1676 NA 16.4411 -0.314 0.75465
## GRADUAT 0.9740 0.6124 0.1626 5.990 2.6e-07 ***
## SF_RATIO -1.5182 -0.2762 0.5621 -2.701 0.00952 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.35 on 48 degrees of freedom
## (2 пропущенных наблюдений удалены)
## Multiple R-squared: 0.6133, Adjusted R-squared: 0.5972
## F-statistic: 38.07 on 2 and 48 DF, p-value: 1.247e-10summary(model_pr_outlier|>lm.beta())##
## Call:
## lm(formula = NEW10 ~ GRADUAT + SF_RATIO, data = pr_data_outlier)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.6995 -8.8741 0.3934 8.0222 31.4411
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -20.0312 NA 14.1534 -1.415 0.163576
## GRADUAT 1.1921 0.6758 0.1443 8.263 1.03e-10 ***
## SF_RATIO -1.9009 -0.3250 0.4784 -3.973 0.000242 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.87 on 47 degrees of freedom
## (2 пропущенных наблюдений удалены)
## Multiple R-squared: 0.7332, Adjusted R-squared: 0.7218
## F-statistic: 64.57 on 2 and 47 DF, p-value: 3.284e-14Комментарий: замечаем существенную разницу в R^2, p-value и значимости SF_RATIO.
plot(model_pr_outlier, 1:5)
AIC(model_pr_full)## [1] 428.221AIC(model_pr_outlier)## [1] 402.2787Можем заметить, что модель без выброса лучше по AIC, чем с выбросом.
4.7.10 Итог: прогнозирование
set.seed(456) # Другое зерно для генерации новых значений
max_GRADUAT<-max(pr_data$GRADUAT, na.rm=TRUE)
max_SF_RATIO<-max(pr_data$SF_RATIO, na.rm=TRUE)
new_data <- data.frame(
GRADUAT = max_GRADUAT - 1:8,
SF_RATIO = max_SF_RATIO - 1:8
)
# Прогнозируем значение NEW10
predicted_values <- predict(model_pr_outlier, newdata = new_data)
cat("Прогнозируемые значения NEW10:\n")## Прогнозируемые значения NEW10:print(predicted_values)## 1 2 3 4 5 6 7 8
## 59.72933 60.43809 61.14684 61.85560 62.56435 63.27311 63.98186 64.69061# 6. Доверительные и предсказательные интервалы
# Доверительный интервал (confidence interval) для среднего значения NEW10
confidence_intervals <- predict(model_pr_outlier, newdata = new_data, interval = "confidence", level = 0.95)
print("Доверительные интервалы (95%):\n")## [1] "Доверительные интервалы (95%):\n"print(confidence_intervals)## fit lwr upr
## 1 59.72933 46.77615 72.68252
## 2 60.43809 48.53109 72.34508
## 3 61.14684 50.27377 72.01991
## 4 61.85560 52.00033 71.71086
## 5 62.56435 53.70520 71.42350
## 6 63.27311 55.38017 71.16604
## 7 63.98186 57.01279 70.95093
## 8 64.69061 58.58381 70.79742# Предсказательный интервал (prediction interval) для отдельного наблюдения NEW10
prediction_intervals <- predict(model_pr_outlier, newdata = new_data, interval = "prediction", level = 0.95)
print("Предсказательные интервалы (95%):\n")## [1] "Предсказательные интервалы (95%):\n"print(prediction_intervals)## fit lwr upr
## 1 59.72933 30.78140 88.67727
## 2 60.43809 31.94292 88.93325
## 3 61.14684 33.06800 89.22569
## 4 61.85560 34.15498 89.55621
## 5 62.56435 35.20229 89.92641
## 6 63.27311 36.20843 90.33778
## 7 63.98186 37.17205 90.79167
## 8 64.69061 38.09191 91.28932# Создаем data.frame для прогнозов и интервалов
predictions_df <- data.frame(
GRADUAT = new_data$GRADUAT,
SF_RATIO = new_data$SF_RATIO,
Predicted = predicted_values,
Lower_Conf = confidence_intervals[, "lwr"],
Upper_Conf = confidence_intervals[, "upr"],
Lower_Pred = prediction_intervals[, "lwr"],
Upper_Pred = prediction_intervals[, "upr"]
)
predictions_df# Визуализация (только для демонстрации, нужен более сложный график для нескольких предикторов)
ggplot(predictions_df, aes(x = 1:nrow(predictions_df), y = Predicted)) + # Условная ось x, т.к. несколько предикторов
geom_point() +
geom_errorbar(aes(ymin = Lower_Conf, ymax = Upper_Conf), color = "blue", width = 0.2) + # Доверительные интервалы
geom_errorbar(aes(ymin = Lower_Pred, ymax = Upper_Pred), color = "red", width = 0.2) + # Предсказательные интервалы
labs(title = "Прогнозы NEW10 с доверительными и предсказательными интервалами",
x = "Прогноз (номер)",
y = "NEW10") +
scale_x_continuous(breaks = 1:nrow(predictions_df)) + # Подписи для оси x
theme_bw()
# annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.9, label = "Синие: Доверительные интервалы", color = "blue", hjust = 0) +
# annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.8, label = "Красные: Предсказательные интервалы", color = "red", hjust = 0)4.8 Анализ общественных вузов
pub_data<-col_I_regr_split$public
pub_data<-pub_data|>mutate(NEW10 = log(NEW10))|>rename(log_NEW10 = NEW10)
print(paste0("Число наблюдений: ", dim(pub_data)[1]))## [1] "Число наблюдений: 123"colSums(is.na(pub_data))## ...1 log_ADD_FEE BOOK R_B_COST OUT_STAT log_NEW10
## 0 33 2 0 1 15
## SF_RATIO PH_D GRADUAT PPIND
## 0 5 4 0Много пропусков по log_ADD_FEE (около четверти).
4.8.1 График
plot_numeric_pairs(pub_data, "Общественные вузы")
4.8.2 Выбросы “на глаз”
outlier_pub_out_eye<-pub_data[pub_data$OUT_STAT>14000,]
outlier_pub_out_eyeplot_numeric_pairs(pub_data, "Общественные вузы", list(red = outlier_pub_out_eye))
4.8.3 Линейная модель
library(lm.beta)
model_pub_full <- lm(log_NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT,
data = pub_data)
summary(model_pub_full|>lm.beta())##
## Call:
## lm(formula = log_NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D +
## SF_RATIO + GRADUAT + OUT_STAT, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.23717 -0.22303 0.04705 0.22529 0.89431
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -4.821e-02 NA 8.271e-01 -0.058 0.9537
## BOOK 8.739e-04 1.632e-01 5.489e-04 1.592 0.1161
## log_ADD_FEE 4.078e-03 7.717e-03 5.368e-02 0.076 0.9397
## R_B_COST 8.749e-05 1.289e-01 7.275e-05 1.203 0.2334
## PH_D 2.303e-02 2.706e-01 8.974e-03 2.566 0.0126 *
## SF_RATIO -9.484e-03 -7.112e-02 1.289e-02 -0.736 0.4645
## GRADUAT 1.970e-02 5.464e-01 4.132e-03 4.768 1.07e-05 ***
## OUT_STAT -5.986e-05 -2.567e-01 2.867e-05 -2.088 0.0407 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4196 on 66 degrees of freedom
## (49 пропущенных наблюдений удалены)
## Multiple R-squared: 0.4814, Adjusted R-squared: 0.4264
## F-statistic: 8.752 on 7 and 66 DF, p-value: 1.449e-074.8.4 Доверительные эллипсы
plot_confidence_ellipse(model_pub_full, c(2, 3))
plot_confidence_ellipse(model_pub_full, c(4, 5))
plot_confidence_ellipse(model_pub_full, c(6, 7))
plot_confidence_ellipse(model_pub_full, c(8, 2))
# Автор: -Евгений-
# Вычисление множественной корреляции
calc_multicollinearity <- function(var, data) {
X <- as.matrix(data[, setdiff(names(data), var)]) # Матрица остальных регрессоров
y <- as.matrix(data[[var]]) # Текущий регрессор как зависимая переменная
model <- lm(y ~ X)
sqrt(summary(model)$r.squared)
}
calc_partial_correlation <- function(var, data) {
# Остатки регрессии текущего регрессора на остальные регрессоры
data<-data|>na.omit()
model_X <- lm(data[[var]] ~ ., data = data[, !names(data) %in% c(var, "log_NEW10")])
res_X <- residuals(model_X)
# Остатки регрессии Y на остальные регрессоры
model_Y <- lm(log_NEW10 ~ ., data = data[, !names(data) %in% var])
res_Y <- residuals(model_Y)
cor(res_X, res_Y, use = "complete.obs") # Корреляция остатков
}
regressors <- setdiff(names(pub_data), c("...1", "PPIND", "log_NEW10"))
regressors_full<-c(regressors, "log_NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pub_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pub_data[,regressors_full])
)
print(redundancy_table)## Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE 0.4883582 0.009351778
## BOOK BOOK 0.5291391 0.192313348
## R_B_COST R_B_COST 0.5747925 0.146434350
## OUT_STAT OUT_STAT 0.7158033 -0.248923976
## SF_RATIO SF_RATIO 0.4073614 -0.090198497
## PH_D PH_D 0.5981262 0.301176912
## GRADUAT GRADUAT 0.7449868 0.506130402Исключим log_ADD_FEE и SF_RATIO (так как у них низкие частные корреляции)
AIC(model_pub_full)## [1] 90.99499summary(model_pub_full)##
## Call:
## lm(formula = log_NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D +
## SF_RATIO + GRADUAT + OUT_STAT, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.23717 -0.22303 0.04705 0.22529 0.89431
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.821e-02 8.271e-01 -0.058 0.9537
## BOOK 8.739e-04 5.489e-04 1.592 0.1161
## log_ADD_FEE 4.078e-03 5.368e-02 0.076 0.9397
## R_B_COST 8.749e-05 7.275e-05 1.203 0.2334
## PH_D 2.303e-02 8.974e-03 2.566 0.0126 *
## SF_RATIO -9.484e-03 1.289e-02 -0.736 0.4645
## GRADUAT 1.970e-02 4.132e-03 4.768 1.07e-05 ***
## OUT_STAT -5.986e-05 2.867e-05 -2.088 0.0407 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4196 on 66 degrees of freedom
## (49 пропущенных наблюдений удалены)
## Multiple R-squared: 0.4814, Adjusted R-squared: 0.4264
## F-statistic: 8.752 on 7 and 66 DF, p-value: 1.449e-07model_pub_nono<-lm(log_NEW10 ~ R_B_COST + PH_D + GRADUAT + OUT_STAT + BOOK, data=pub_data)
summary(model_pub_full|>lm.beta())##
## Call:
## lm(formula = log_NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D +
## SF_RATIO + GRADUAT + OUT_STAT, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.23717 -0.22303 0.04705 0.22529 0.89431
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -4.821e-02 NA 8.271e-01 -0.058 0.9537
## BOOK 8.739e-04 1.632e-01 5.489e-04 1.592 0.1161
## log_ADD_FEE 4.078e-03 7.717e-03 5.368e-02 0.076 0.9397
## R_B_COST 8.749e-05 1.289e-01 7.275e-05 1.203 0.2334
## PH_D 2.303e-02 2.706e-01 8.974e-03 2.566 0.0126 *
## SF_RATIO -9.484e-03 -7.112e-02 1.289e-02 -0.736 0.4645
## GRADUAT 1.970e-02 5.464e-01 4.132e-03 4.768 1.07e-05 ***
## OUT_STAT -5.986e-05 -2.567e-01 2.867e-05 -2.088 0.0407 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4196 on 66 degrees of freedom
## (49 пропущенных наблюдений удалены)
## Multiple R-squared: 0.4814, Adjusted R-squared: 0.4264
## F-statistic: 8.752 on 7 and 66 DF, p-value: 1.449e-07AIC(model_pub_full)## [1] 90.99499summary(model_pub_nono|>lm.beta())##
## Call:
## lm(formula = log_NEW10 ~ R_B_COST + PH_D + GRADUAT + OUT_STAT +
## BOOK, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.21996 -0.24955 -0.00251 0.24842 0.81900
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 7.390e-02 NA 5.501e-01 0.134 0.8934
## R_B_COST 1.272e-04 1.806e-01 6.475e-05 1.964 0.0525 .
## PH_D 1.748e-02 2.156e-01 7.462e-03 2.342 0.0213 *
## GRADUAT 1.808e-02 4.905e-01 3.432e-03 5.268 9.03e-07 ***
## OUT_STAT -5.027e-05 -2.003e-01 2.389e-05 -2.104 0.0381 *
## BOOK 1.044e-03 1.903e-01 4.510e-04 2.315 0.0228 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4225 on 92 degrees of freedom
## (25 пропущенных наблюдений удалены)
## Multiple R-squared: 0.4772, Adjusted R-squared: 0.4488
## F-statistic: 16.8 on 5 and 92 DF, p-value: 9.211e-12AIC(model_pub_nono)## [1] 117.0556Комментарий: AIC увеличился, p-value уменьшился. Теперь каждый признак, кроме R_B_COST является значимым (при alpha=0.05).
regressors <- setdiff(names(pub_data), c("...1", "PPIND", "log_NEW10", "log_ADD_FEE", "SF_RATIO"))
regressors_full<-c(regressors, "log_NEW10")
redundancy_table <- data.frame(
Регрессор = regressors,
Множественная_корреляция = sapply(regressors, calc_multicollinearity, pub_data[, regressors_full]),
Частная_корреляция = sapply(regressors, calc_partial_correlation, pub_data[,regressors_full])
)
print(redundancy_table)## Регрессор Множественная_корреляция Частная_корреляция
## BOOK BOOK 0.4535409 0.2346185
## R_B_COST R_B_COST 0.5958627 0.2006126
## OUT_STAT OUT_STAT 0.6338016 -0.2142731
## PH_D PH_D 0.6059169 0.2372086
## GRADUAT GRADUAT 0.7046461 0.4813910Тут уже не так понятно. Перейдём к пошаговому AIC.
AIC(model_pub_full)## [1] 90.99499AIC(model_pub_nono)## [1] 117.05564.8.5 OLS. Пошаговый AIC.
4.8.5.1 Backward
library(olsrr)
backward_aic<-ols_step_backward_aic(model_pub_full, progress=TRUE, details=TRUE)## Backward Elimination Method
## ---------------------------
##
## Candidate Terms:
##
## 1. BOOK
## 2. log_ADD_FEE
## 3. R_B_COST
## 4. PH_D
## 5. SF_RATIO
## 6. GRADUAT
## 7. OUT_STAT
##
##
## Step => 0
## Model => log_NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 90.99499
##
## Initiating stepwise selection...
##
## Table: Removing Existing Variables
## -------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## log_ADD_FEE 1 89.001 107.434 -119.335 0.48135 0.43490
## SF_RATIO 1 89.599 108.032 -118.859 0.47714 0.43032
## R_B_COST 1 90.599 109.032 -118.062 0.47003 0.42257
## BOOK 1 91.784 110.216 -117.115 0.46148 0.41325
## OUT_STAT 1 93.728 112.161 -115.555 0.44714 0.39763
## PH_D 1 96.032 114.464 -113.699 0.42966 0.37858
## GRADUAT 1 110.895 129.327 -101.488 0.30279 0.24035
## -------------------------------------------------------------------------
##
## Step => 1
## Removed => log_ADD_FEE
## Model => log_NEW10 ~ BOOK + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT
## AIC => 89.00146
##
## Table: Removing Existing Variables
## -----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -----------------------------------------------------------------------
## SF_RATIO 1 87.604 103.733 -121.080 0.47711 0.43866
## R_B_COST 1 88.775 104.904 -120.111 0.46877 0.42970
## BOOK 1 89.835 105.964 -119.232 0.46110 0.42148
## OUT_STAT 1 91.997 108.125 -117.436 0.44513 0.40433
## PH_D 1 94.185 110.314 -115.609 0.42847 0.38645
## GRADUAT 1 111.437 127.565 -100.972 0.27842 0.22537
## -----------------------------------------------------------------------
##
## Step => 2
## Removed => SF_RATIO
## Model => log_NEW10 ~ BOOK + R_B_COST + PH_D + GRADUAT + OUT_STAT
## AIC => 87.60406
##
## Table: Removing Existing Variables
## -----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -----------------------------------------------------------------------
## R_B_COST 1 87.268 101.092 -121.878 0.46522 0.43422
## BOOK 1 88.652 102.476 -120.689 0.45512 0.42354
## OUT_STAT 1 90.007 103.831 -119.523 0.44505 0.41288
## PH_D 1 94.291 108.115 -115.824 0.41198 0.37789
## GRADUAT 1 109.486 123.311 -102.540 0.27794 0.23608
## -----------------------------------------------------------------------
##
## Step => 3
## Removed => R_B_COST
## Model => log_NEW10 ~ BOOK + PH_D + GRADUAT + OUT_STAT
## AIC => 87.26753
##
## Table: Removing Existing Variables
## -----------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -----------------------------------------------------------------------
## OUT_STAT 1 88.405 99.925 -121.264 0.44206 0.41815
## BOOK 1 89.065 100.585 -120.677 0.43706 0.41294
## PH_D 1 96.921 108.442 -113.661 0.37401 0.34718
## GRADUAT 1 109.030 120.550 -102.753 0.26272 0.23112
## -----------------------------------------------------------------------
##
##
## No more variables to be removed.
##
## Variables Removed:
##
## => log_ADD_FEE
## => SF_RATIO
## => R_B_COSTbackward_aic##
##
## Stepwise Summary
## --------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## --------------------------------------------------------------------------
## 0 Full Model 90.995 111.732 -117.098 0.48139 0.42639
## 1 log_ADD_FEE 89.001 107.434 -119.561 0.48135 0.43490
## 2 SF_RATIO 87.604 103.733 -121.356 0.47711 0.43866
## 3 R_B_COST 87.268 101.092 -122.021 0.46522 0.43422
## --------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ---------------------------------------------------------------
## R 0.682 RMSE 0.402
## R-Squared 0.465 MSE 0.162
## Adj. R-Squared 0.434 Coef. Var 12.482
## Pred R-Squared 0.385 AIC 87.268
## MAE 0.300 SBC 101.092
## ---------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 10.423 4 2.606 15.006 0.0000
## Residual 11.981 69 0.174
## Total 22.403 73
## -------------------------------------------------------------------
##
## Parameter Estimates
## ---------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ---------------------------------------------------------------------------------------
## (Intercept) -0.435 0.670 -0.649 0.519 -1.771 0.902
## BOOK 0.001 0.001 0.188 1.906 0.061 0.000 0.002
## PH_D 0.028 0.008 0.329 3.431 0.001 0.012 0.044
## GRADUAT 0.019 0.004 0.538 5.112 0.000 0.012 0.027
## OUT_STAT 0.000 0.000 -0.187 -1.729 0.088 0.000 0.000
## ---------------------------------------------------------------------------------------4.8.5.2 Forward
forward_aic<-ols_step_forward_aic(model_pub_full, progress=TRUE, details=TRUE)## Forward Selection Method
## ------------------------
##
## Candidate Terms:
##
## 1. BOOK
## 2. log_ADD_FEE
## 3. R_B_COST
## 4. PH_D
## 5. SF_RATIO
## 6. GRADUAT
## 7. OUT_STAT
##
##
## Step => 0
## Model => log_NEW10 ~ 1
## AIC => 125.5841
##
## Initiating stepwise selection...
##
## Table: Adding New Variables
## --------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## --------------------------------------------------------------------------
## GRADUAT 1 107.184 114.097 -103.864 0.24094 0.23040
## PH_D 1 107.752 114.664 -103.325 0.23509 0.22447
## R_B_COST 1 119.631 126.543 -92.031 0.10191 0.08943
## BOOK 1 121.310 128.222 -90.431 0.08130 0.06854
## log_ADD_FEE 1 122.064 128.977 -89.711 0.07188 0.05899
## SF_RATIO 1 126.507 133.419 -85.472 0.01445 0.00076
## OUT_STAT 1 127.284 134.197 -84.730 0.00404 -0.00979
## --------------------------------------------------------------------------
##
## Step => 1
## Added => GRADUAT
## Model => log_NEW10 ~ GRADUAT
## AIC => 107.1844
##
## Table: Adding New Variables
## -------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------
## PH_D 1 92.767 101.983 -117.502 0.39197 0.37484
## BOOK 1 96.918 106.134 -113.678 0.35689 0.33877
## OUT_STAT 1 104.070 113.286 -107.072 0.29163 0.27168
## R_B_COST 1 105.322 114.539 -105.913 0.27954 0.25924
## SF_RATIO 1 107.925 117.141 -103.501 0.25375 0.23273
## log_ADD_FEE 1 108.311 117.527 -103.142 0.24984 0.22871
## -------------------------------------------------------------------------
##
## Step => 2
## Added => PH_D
## Model => log_NEW10 ~ GRADUAT + PH_D
## AIC => 92.76671
##
## Table: Adding New Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## BOOK 1 88.405 99.925 -121.264 0.44206 0.41815
## OUT_STAT 1 89.065 100.585 -120.677 0.43706 0.41294
## R_B_COST 1 94.202 105.723 -116.094 0.39659 0.37073
## log_ADD_FEE 1 94.654 106.174 -115.690 0.39290 0.36688
## SF_RATIO 1 94.759 106.280 -115.596 0.39203 0.36597
## ------------------------------------------------------------------------
##
## Step => 3
## Added => BOOK
## Model => log_NEW10 ~ GRADUAT + PH_D + BOOK
## AIC => 88.40477
##
## Table: Adding New Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## OUT_STAT 1 87.268 101.092 -121.878 0.46522 0.43422
## R_B_COST 1 90.007 103.831 -119.523 0.44505 0.41288
## log_ADD_FEE 1 90.016 103.840 -119.515 0.44498 0.41281
## SF_RATIO 1 90.387 104.212 -119.195 0.44219 0.40986
## ------------------------------------------------------------------------
##
## Step => 4
## Added => OUT_STAT
## Model => log_NEW10 ~ GRADUAT + PH_D + BOOK + OUT_STAT
## AIC => 87.26753
##
## Table: Adding New Variables
## ------------------------------------------------------------------------
## Predictor DF AIC SBC SBIC R2 Adj. R2
## ------------------------------------------------------------------------
## R_B_COST 1 87.604 103.733 -121.080 0.47711 0.43866
## SF_RATIO 1 88.775 104.904 -120.111 0.46877 0.42970
## log_ADD_FEE 1 89.109 105.237 -119.834 0.46636 0.42713
## ------------------------------------------------------------------------
##
##
## No more variables to be added.
##
## Variables Selected:
##
## => GRADUAT
## => PH_D
## => BOOK
## => OUT_STATforward_aic##
##
## Stepwise Summary
## --------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## --------------------------------------------------------------------------
## 0 Base Model 125.584 130.192 -85.606 0.00000 0.00000
## 1 GRADUAT 107.184 114.097 -103.864 0.24094 0.23040
## 2 PH_D 92.767 101.983 -117.502 0.39197 0.37484
## 3 BOOK 88.405 99.925 -121.264 0.44206 0.41815
## 4 OUT_STAT 87.268 101.092 -121.878 0.46522 0.43422
## --------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## ---------------------------------------------------------------
## R 0.682 RMSE 0.402
## R-Squared 0.465 MSE 0.162
## Adj. R-Squared 0.434 Coef. Var 12.482
## Pred R-Squared 0.385 AIC 87.268
## MAE 0.300 SBC 101.092
## ---------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 10.423 4 2.606 15.006 0.0000
## Residual 11.981 69 0.174
## Total 22.403 73
## -------------------------------------------------------------------
##
## Parameter Estimates
## ---------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ---------------------------------------------------------------------------------------
## (Intercept) -0.435 0.670 -0.649 0.519 -1.771 0.902
## GRADUAT 0.019 0.004 0.538 5.112 0.000 0.012 0.027
## PH_D 0.028 0.008 0.329 3.431 0.001 0.012 0.044
## BOOK 0.001 0.001 0.188 1.906 0.061 0.000 0.002
## OUT_STAT 0.000 0.000 -0.187 -1.729 0.088 0.000 0.000
## ---------------------------------------------------------------------------------------Комментарий: строим Backward и Forward AIC. Backward AIC исключил всё, кроме SF_RATIO и GRADUAT. Forward AIC включил только GRADUAT и SF_RATIO. Оба пошаговых алгоритма дали один и тот же результат.
Интерпретация параметров: число лучших 10% студентов NEW10 увеличивается с увеличением выпускников GRADUAT, числом докторов (PH_D), ценой за бронирование (BOOK) и уменьшается со стоимостью оплаты OUT_STAT.
Здесь можно сравнить по R^2, например.
4.8.6 Построение и анализ графиков
model_pub_full<-lm(log_NEW10 ~ GRADUAT + PH_D + BOOK + OUT_STAT, data=pub_data)
summary(model_pub_full)##
## Call:
## lm(formula = log_NEW10 ~ GRADUAT + PH_D + BOOK + OUT_STAT, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.21130 -0.24350 -0.01198 0.24097 0.92728
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.967e-02 5.542e-01 -0.108 0.91449
## GRADUAT 1.825e-02 3.483e-03 5.241 9.95e-07 ***
## PH_D 2.271e-02 7.077e-03 3.209 0.00183 **
## BOOK 1.152e-03 4.545e-04 2.535 0.01290 *
## OUT_STAT -3.446e-05 2.284e-05 -1.509 0.13474
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4289 on 93 degrees of freedom
## (25 пропущенных наблюдений удалены)
## Multiple R-squared: 0.4553, Adjusted R-squared: 0.4319
## F-statistic: 19.43 on 4 and 93 DF, p-value: 1.196e-11AIC(model_pub_full)## [1] 119.0812summary(model_pub_nono)##
## Call:
## lm(formula = log_NEW10 ~ R_B_COST + PH_D + GRADUAT + OUT_STAT +
## BOOK, data = pub_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.21996 -0.24955 -0.00251 0.24842 0.81900
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.390e-02 5.501e-01 0.134 0.8934
## R_B_COST 1.272e-04 6.475e-05 1.964 0.0525 .
## PH_D 1.748e-02 7.462e-03 2.342 0.0213 *
## GRADUAT 1.808e-02 3.432e-03 5.268 9.03e-07 ***
## OUT_STAT -5.027e-05 2.389e-05 -2.104 0.0381 *
## BOOK 1.044e-03 4.510e-04 2.315 0.0228 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4225 on 92 degrees of freedom
## (25 пропущенных наблюдений удалены)
## Multiple R-squared: 0.4772, Adjusted R-squared: 0.4488
## F-statistic: 16.8 on 5 and 92 DF, p-value: 9.211e-12AIC(model_pub_nono)## [1] 117.0556Комментарий: эта модель в сравнении с model_pub_nono имеет более высокий p-value и более низкий R^2.
plot(model_pub_full, 1:5)
4.8.7 Поиск outliers по Residuals vs. Deleted Residuals
plot(residuals(model_pub_full), rstudent(model_pub_full))
Комментарий: особых выбросов нет.
influence.measures(model_pub_full)## Influence measures of
## lm(formula = log_NEW10 ~ GRADUAT + PH_D + BOOK + OUT_STAT, data = pub_data) :
##
## dfb.1_ dfb.GRAD dfb.PH_D dfb.BOOK dfb.OUT_ dffit cov.r cook.d
## 1 -0.042131 -0.133457 3.02e-02 0.016099 1.27e-01 -0.17853 1.047 6.39e-03
## 3 0.103254 0.106100 -1.30e-01 -0.039284 4.56e-02 -0.20747 1.147 8.67e-03
## 5 0.025264 0.031849 -1.62e-02 -0.037843 -1.18e-02 -0.07080 1.072 1.01e-03
## 6 0.042930 -0.026515 -3.75e-02 0.020769 2.19e-03 0.07699 1.074 1.20e-03
## 7 -0.001421 -0.001482 2.05e-03 -0.000293 3.65e-05 0.00295 1.083 1.76e-06
## 8 0.256242 -0.005185 -1.43e-01 -0.095065 -1.30e-01 0.30174 1.034 1.81e-02
## 9 -0.223857 0.071406 9.59e-02 0.076188 2.17e-01 0.40619 0.932 3.22e-02
## 11 -0.350149 -0.136018 1.30e-01 0.312988 3.54e-01 0.52495 1.037 5.43e-02
## 12 -0.121942 0.075940 3.86e-02 0.138385 -2.94e-02 0.21075 1.137 8.94e-03
## 13 -0.215844 -0.020160 1.53e-01 0.134773 -1.49e-03 0.27751 1.055 1.54e-02
## 15 -0.321287 -0.077636 2.49e-01 0.025996 2.44e-01 0.44045 0.919 3.77e-02
## 16 -0.250015 -0.049278 2.92e-01 -0.032726 -4.84e-02 0.38214 0.833 2.80e-02
## 17 -0.081306 0.044052 6.07e-02 0.007316 -2.62e-02 -0.11473 1.077 2.65e-03
## 18 0.003194 -0.000841 1.97e-02 -0.041553 -4.86e-03 0.05185 1.085 5.43e-04
## 19 0.009679 0.007517 -7.39e-03 0.003912 -2.18e-02 -0.02766 1.122 1.55e-04
## 20 0.141867 0.108893 -1.32e-01 0.020169 -1.07e-01 0.23226 0.956 1.06e-02
## 21 -0.012381 0.071752 1.14e-02 0.006585 -6.69e-02 0.12529 1.036 3.15e-03
## 22 0.036851 -0.032196 3.27e-02 -0.086552 -3.76e-02 0.12703 1.053 3.24e-03
## 23 -0.112694 0.134024 1.63e-02 0.188038 -1.08e-01 0.30489 1.045 1.85e-02
## 24 0.134601 0.381840 -1.22e-01 -0.241573 -1.15e-01 -0.53441 0.843 5.46e-02
## 25 0.014454 0.059799 3.86e-02 -0.077258 -1.06e-01 0.13967 1.084 3.93e-03
## 26 0.006213 0.014783 -3.94e-03 0.011396 -3.38e-02 0.04477 1.100 4.05e-04
## 27 -0.003111 0.012714 4.84e-02 -0.119863 2.89e-03 -0.13920 1.130 3.91e-03
## 28 -0.224463 -0.046283 1.85e-01 0.043128 1.35e-02 -0.27959 0.924 1.53e-02
## 29 -0.199615 -0.064122 1.55e-01 0.074507 2.84e-02 -0.22725 1.057 1.03e-02
## 31 0.024836 0.118218 9.65e-05 -0.070988 -8.35e-02 0.15237 1.108 4.68e-03
## 32 0.010801 -0.068147 -4.39e-04 0.018210 2.63e-02 0.08799 1.083 1.56e-03
## 33 0.546814 0.131824 -4.82e-01 -0.062290 -6.47e-02 0.58654 0.888 6.62e-02
## 34 -0.017402 -0.014653 2.72e-02 -0.010574 -9.31e-03 -0.03488 1.106 2.46e-04
## 35 0.000434 -0.007105 1.32e-03 0.002973 -6.78e-03 -0.02361 1.072 1.13e-04
## 36 -0.030540 -0.044774 4.61e-02 -0.026909 1.86e-02 -0.07377 1.075 1.10e-03
## 38 0.038315 0.014707 -3.98e-02 0.010623 -9.85e-03 0.05040 1.086 5.13e-04
## 39 0.031172 -0.009854 -1.04e-02 -0.062093 2.45e-02 -0.11232 1.044 2.53e-03
## 41 -0.001574 0.010763 -4.89e-03 0.005513 -4.69e-04 -0.01433 1.110 4.15e-05
## 42 0.033069 -0.087660 1.88e-02 -0.011504 -2.87e-02 0.15290 1.040 4.69e-03
## 44 0.043104 0.028398 -9.65e-02 0.086072 1.48e-02 -0.12728 1.088 3.26e-03
## 45 0.027188 -0.009637 -5.15e-02 0.053262 -1.60e-03 -0.11268 1.036 2.55e-03
## 46 0.032962 -0.006529 -6.83e-03 -0.044942 8.94e-03 0.07515 1.064 1.14e-03
## 47 -0.130253 0.121583 6.48e-02 0.136264 -2.07e-01 -0.38168 0.845 2.80e-02
## 48 -0.000454 -0.146705 -1.18e-01 0.270710 8.78e-02 -0.40760 0.809 3.17e-02
## 49 0.091508 -0.006039 -1.23e-01 0.115218 -6.34e-02 -0.22342 1.046 9.98e-03
## 51 0.020780 0.007590 -8.55e-03 -0.017144 -1.30e-02 0.03116 1.082 1.96e-04
## 52 0.005839 -0.001717 -4.46e-03 -0.003079 4.69e-03 0.01013 1.100 2.07e-05
## 53 -0.002512 -0.001853 1.40e-03 0.000344 3.86e-03 -0.00499 1.094 5.04e-06
## 54 0.017187 0.003092 -5.29e-03 -0.012185 -1.72e-02 0.02580 1.111 1.35e-04
## 55 -0.007263 -0.006925 6.30e-03 0.000973 1.36e-02 0.03071 1.066 1.91e-04
## 56 -0.102366 -0.142662 4.14e-02 0.142251 1.49e-01 0.25798 0.993 1.32e-02
## 58 0.071336 0.007929 -5.81e-02 -0.012835 -6.21e-03 0.07724 1.110 1.20e-03
## 59 -0.020325 0.005678 1.20e-02 0.003531 7.44e-03 -0.02635 1.101 1.40e-04
## 61 0.026528 0.015272 3.82e-02 -0.084232 -8.32e-02 0.12538 1.081 3.17e-03
## 63 0.087394 -0.059968 1.25e-02 -0.090035 -1.38e-01 -0.24911 1.013 1.23e-02
## 64 -0.170093 -0.191554 2.91e-01 -0.103360 -5.94e-02 0.34442 1.146 2.38e-02
## 66 0.247524 -0.418324 -1.60e-01 -0.152779 3.84e-01 -0.69659 0.693 8.92e-02
## 67 0.034292 0.137055 -7.46e-02 0.059278 -8.08e-02 0.16454 1.130 5.46e-03
## 68 0.026783 -0.003626 -5.07e-02 0.024273 4.16e-02 -0.08292 1.068 1.39e-03
## 69 -0.019209 -0.001041 2.31e-02 -0.007225 -4.35e-04 0.03438 1.075 2.39e-04
## 70 0.012013 0.239901 -9.68e-03 -0.078056 -1.17e-01 0.28916 1.023 1.66e-02
## 71 -0.092085 -0.005378 1.14e-01 -0.044926 -4.20e-02 -0.14123 1.070 4.01e-03
## 72 0.093203 0.066536 -2.64e-01 0.242285 1.34e-01 -0.34264 1.152 2.35e-02
## 73 -0.158831 -0.232355 1.86e-01 -0.033745 1.25e-01 -0.32859 0.894 2.10e-02
## 74 -0.110318 0.092402 9.43e-03 0.127219 -2.19e-02 -0.23596 0.963 1.10e-02
## 75 -0.062227 -0.055431 5.68e-02 0.012618 2.93e-02 -0.09321 1.070 1.75e-03
## 77 0.002545 0.003570 -4.48e-03 0.002996 -2.57e-03 -0.00797 1.079 1.28e-05
## 79 0.000397 -0.000864 1.38e-02 -0.043803 1.39e-02 -0.05684 1.119 6.53e-04
## 80 -0.003717 -0.005744 -5.30e-04 0.022151 -1.00e-02 0.03494 1.106 2.47e-04
## 83 0.030763 -0.027421 -4.54e-02 0.016224 7.42e-02 0.10076 1.086 2.05e-03
## 84 0.051012 0.143654 -5.12e-02 -0.013781 -1.43e-01 -0.19210 1.059 7.40e-03
## 85 0.023201 0.012497 -5.37e-02 0.070858 -2.02e-02 -0.09385 1.147 1.78e-03
## 86 -0.045718 0.003991 1.31e-02 0.126471 -1.41e-01 -0.29579 0.921 1.71e-02
## 87 -0.009540 -0.062359 6.97e-02 -0.098075 2.55e-03 -0.12517 1.147 3.16e-03
## 88 0.014200 0.014410 -2.16e-04 -0.027171 -1.11e-02 0.04109 1.078 3.41e-04
## 90 0.008670 -0.012445 1.30e-02 -0.071909 3.83e-02 -0.10831 1.080 2.36e-03
## 91 0.029025 -0.048617 7.56e-02 -0.120907 -9.13e-02 0.17578 1.144 6.23e-03
## 92 0.430023 -0.528612 -4.13e-01 0.206695 3.71e-01 0.83959 1.084 1.38e-01
## 93 0.045189 -0.017476 1.07e-02 -0.050952 -6.48e-02 0.10855 1.083 2.37e-03
## 96 0.196067 -0.069055 -8.44e-02 -0.082115 -9.60e-02 0.26129 1.081 1.37e-02
## 97 -0.009366 0.027867 1.62e-02 -0.002606 -4.52e-02 0.05736 1.107 6.65e-04
## 98 0.027003 0.036114 -7.06e-03 -0.057481 -1.77e-02 -0.07484 1.166 1.13e-03
## 100 0.075991 0.026069 1.59e-02 0.035124 -3.98e-01 -0.48472 1.107 4.67e-02
## 101 -0.028942 0.123097 -2.87e-02 0.017744 3.96e-02 0.18388 1.111 6.80e-03
## 102 -0.004197 0.007296 5.48e-03 -0.003717 -9.50e-03 -0.01427 1.096 4.12e-05
## 103 0.007016 -0.020480 -1.38e-03 -0.009095 2.28e-02 0.03411 1.105 2.35e-04
## 104 0.123240 0.179060 -1.41e-01 0.026554 -1.63e-01 -0.24651 1.098 1.22e-02
## 105 0.023884 -0.063654 8.23e-02 -0.146901 -7.63e-02 -0.20001 1.083 8.03e-03
## 106 0.055488 -0.001506 -4.34e-02 -0.022109 5.94e-04 -0.07177 1.090 1.04e-03
## 107 0.004810 0.000641 1.61e-02 -0.041155 -1.16e-02 -0.04542 1.121 4.17e-04
## 109 0.015535 -0.011063 -2.23e-02 0.019857 3.92e-03 -0.03938 1.093 3.13e-04
## 110 -0.045418 0.099942 3.16e-02 -0.007424 -8.10e-02 -0.14449 1.071 4.20e-03
## 111 0.003248 0.005327 -1.45e-03 -0.004667 -4.85e-03 -0.00846 1.105 1.45e-05
## 113 0.037474 -0.129469 3.64e-03 -0.054592 7.57e-02 -0.17693 1.065 6.28e-03
## 114 0.223708 0.055697 -2.56e-01 -0.035107 2.07e-01 0.42656 0.901 3.53e-02
## 115 -0.039680 0.195298 2.60e-03 -0.087338 5.77e-02 0.32221 1.088 2.07e-02
## 116 -0.012937 -0.026990 1.32e-02 0.008373 9.61e-03 -0.03206 1.140 2.08e-04
## 117 -0.137025 0.039881 4.03e-02 -0.106701 4.40e-01 0.61275 1.075 7.40e-02
## 118 0.001980 -0.000154 -1.45e-03 -0.002237 1.56e-03 -0.00497 1.084 4.99e-06
## 119 -0.195111 0.041427 1.92e-01 -0.064464 -1.77e-02 -0.25500 1.043 1.30e-02
## 122 -0.049622 -0.115513 9.73e-02 -0.034789 2.04e-02 0.14702 1.104 4.35e-03
## 123 0.002463 0.058149 1.39e-02 -0.020111 -7.73e-02 0.09267 1.112 1.73e-03
## hat inf
## 1 0.0367
## 3 0.0988
## 5 0.0252
## 6 0.0278
## 7 0.0254
## 8 0.0578
## 9 0.0482
## 11 0.1039
## 12 0.0931
## 13 0.0613
## 15 0.0513
## 16 0.0283 *
## 17 0.0378
## 18 0.0312
## 19 0.0595
## 20 0.0232
## 21 0.0210
## 22 0.0278
## 23 0.0626
## 24 0.0515
## 25 0.0471
## 26 0.0428
## 27 0.0772
## 28 0.0256
## 29 0.0513
## 31 0.0647
## 32 0.0365
## 33 0.0693
## 34 0.0470
## 35 0.0171
## 36 0.0278
## 38 0.0325
## 39 0.0210
## 41 0.0491
## 42 0.0283
## 44 0.0472
## 45 0.0181
## 46 0.0214
## 47 0.0296
## 48 0.0288 *
## 49 0.0460
## 51 0.0267
## 52 0.0403
## 53 0.0354
## 54 0.0508
## 55 0.0137
## 56 0.0355
## 58 0.0543
## 59 0.0418
## 61 0.0420
## 63 0.0393
## 64 0.1210
## 66 0.0505 *
## 67 0.0814
## 68 0.0257
## 69 0.0209
## 70 0.0510
## 71 0.0391
## 72 0.1239
## 73 0.0286
## 74 0.0248
## 75 0.0290
## 77 0.0218
## 79 0.0593
## 80 0.0472
## 83 0.0403
## 84 0.0447
## 85 0.0845
## 86 0.0277
## 87 0.0877
## 88 0.0242
## 90 0.0379
## 91 0.0922
## 92 0.1798 *
## 93 0.0399
## 96 0.0703
## 97 0.0495
## 98 0.0974 *
## 100 0.1264
## 101 0.0723
## 102 0.0374
## 103 0.0464
## 104 0.0760
## 105 0.0584
## 106 0.0386
## 107 0.0597
## 109 0.0362
## 110 0.0403
## 111 0.0452
## 113 0.0442
## 114 0.0452
## 115 0.0860
## 116 0.0745
## 117 0.1361
## 118 0.0261
## 119 0.0515
## 122 0.0610
## 123 0.05794.8.8 Поиск outliers по Cook и Mahalanobis
4.8.8.1 Cook
# Рассчитаем Cook's distance для модели с логарифмами
cooks_vals <- cooks.distance(model_pub_full)
cook_df <- data.frame(
obs = row.names(model_pub_full$model),
cook = cooks_vals
) %>%
arrange(desc(cook))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(cook_df))
plot(
x, cook_df$cook,
type = "h",
main = "Cook's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Cook's distance",
xaxt = "n" # убираем автоматическую ось X
)
axis(side = 1, at = x, labels = cook_df$obs, las = 2, cex.axis = 0.8)
abline(h = 1, col = "red", lty = 2)
Комментарий: высокие значения у 29 и 38.
4.8.8.2 Mahalanobis
# 7.2. Расстояние Махаланобиса (Mahalanobis Distance)
# Расстояние Махаланобиса в пространстве предикторов (регрессоров)
mahalanobis_distance <- mahalanobis(model.matrix(model_pub_full)[,-1],
center = colMeans(model.matrix(model_pub_full)[,-1]),
cov = cov(model.matrix(model_pub_full)[,-1]))
mah_df <- data.frame(
obs = row.names(model_pub_full$model),
mah = mahalanobis_distance
) %>%
arrange(desc(mah))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(mah_df))
alpha<-0.05
# Scatter plot для расстояния Махаланобиса
plot(
x, mah_df$mah,
type = "h",
main = "Mahalanobis's Distance",
xlab = "Номер наблюдения (из исходного df)",
ylab = "Mahalanobis's distance",
xaxt = "n"
)
axis(side = 1, at = x, labels = mah_df$obs, las = 2, cex.axis = 0.8)
abline(h = qchisq(p = 1 - alpha, df = nrow(mah_df)), col = "red", lty = 2)
4.8.9 Итог: прогнозирование
set.seed(456) # Другое зерно для генерации новых значений
max_GRADUAT<-max(pub_data$GRADUAT, na.rm=TRUE)
max_BOOK<-max(pub_data$BOOK, na.rm=TRUE)
max_OUT_STAT<-max(pub_data$OUT_STAT, na.rm=TRUE)
max_PH_D<-max(pub_data$PH_D, na.rm=TRUE)
new_data <- data.frame(
GRADUAT = max_GRADUAT - 1:8,
BOOK = max_BOOK - 1:8,
OUT_STAT = max_OUT_STAT - seq(100, 800, by=100),
PH_D = max_PH_D - 1:8
)
# Прогнозируем значение NEW10
predicted_values <- predict(model_pub_full, newdata = new_data)
cat("Прогнозируемые значения NEW10:\n")## Прогнозируемые значения NEW10:print(predicted_values)## 1 2 3 4 5 6 7 8
## 4.330459 4.291791 4.253122 4.214454 4.175785 4.137117 4.098448 4.059780# 6. Доверительные и предсказательные интервалы
# Доверительный интервал (confidence interval) для среднего значения NEW10
confidence_intervals <- predict(model_pub_full, newdata = new_data, interval = "confidence", level = 0.95)
print("Доверительные интервалы (95%):\n")## [1] "Доверительные интервалы (95%):\n"print(confidence_intervals)## fit lwr upr
## 1 4.330459 3.917940 4.742979
## 2 4.291791 3.883437 4.700145
## 3 4.253122 3.848577 4.657668
## 4 4.214454 3.813349 4.615558
## 5 4.175785 3.777744 4.573826
## 6 4.137117 3.741753 4.532480
## 7 4.098448 3.705369 4.491527
## 8 4.059780 3.668584 4.450975# Предсказательный интервал (prediction interval) для отдельного наблюдения NEW10
prediction_intervals <- predict(model_pub_full, newdata = new_data, interval = "prediction", level = 0.95)
print("Предсказательные интервалы (95%):\n")## [1] "Предсказательные интервалы (95%):\n"print(prediction_intervals)## fit lwr upr
## 1 4.330459 3.384023 5.276896
## 2 4.291791 3.347162 5.236419
## 3 4.253122 3.310134 5.196111
## 4 4.214454 3.272936 5.155971
## 5 4.175785 3.235569 5.116001
## 6 4.137117 3.198031 5.076202
## 7 4.098448 3.160322 5.036574
## 8 4.059780 3.122441 4.997118# Создаем data.frame для прогнозов и интервалов
predictions_df <- data.frame(
GRADUAT = new_data$GRADUAT,
OUT_STAT = new_data$OUT_STAT,
BOOK = new_data$BOOK,
PH_D = new_data$PH_D,
Predicted = predicted_values,
Lower_Conf = confidence_intervals[, "lwr"],
Upper_Conf = confidence_intervals[, "upr"],
Lower_Pred = prediction_intervals[, "lwr"],
Upper_Pred = prediction_intervals[, "upr"]
)
predictions_df# Визуализация (только для демонстрации, нужен более сложный график для нескольких предикторов)
ggplot(predictions_df, aes(x = 1:nrow(predictions_df), y = Predicted)) + # Условная ось x, т.к. несколько предикторов
geom_point() +
geom_errorbar(aes(ymin = Lower_Conf, ymax = Upper_Conf), color = "blue", width = 0.2) + # Доверительные интервалы
geom_errorbar(aes(ymin = Lower_Pred, ymax = Upper_Pred), color = "red", width = 0.2) + # Предсказательные интервалы
labs(title = "Прогнозы NEW10 с доверительными и предсказательными интервалами",
x = "Прогноз (номер)",
y = "NEW10") +
scale_x_continuous(breaks = 1:nrow(predictions_df)) + # Подписи для оси x
theme_bw()
# annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.9, label = "Синие: Доверительные интервалы", color = "blue", hjust = 0) +
# annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.8, label = "Красные: Предсказательные интервалы", color = "red", hjust = 0)